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FORWARD 


The attainable service life of a turbine blade is very sensitive to the amplitude of vibratory 
stress which, under resonant conditions, is inversely proportional to the amount of 
damping in the system. Inherent damping in turbine-bladed disk assemblies is relatively 
small. Therefore, the introduction of effective damping is instrumental in reducing 
alternating stress and extending blade life. 

The most commonly-used technique to increase the damping in turbine blades is to provide 
for vibrational energy dissipation, through the use of friction interfaces between adjacent 
blades. Many different damping configurations, based on this concept, have been used in 
the past; the designs of which were largely based on empirical data and on experience with 
previous applications. There have been several published attempts to analytically evaluate 
these designs. For the most part, they consist of simple studies, aimed at understanding 
the damping phenomenon, and do not provide useful design tools. Therefore, the 
objectives of this program are to develop and to experimentally validate the analytical tools 
necessary for the evaluation of friction dampers on a detail design level. Secondly, it is 
hoped that the extensive experimental database generated in this program can be used by 
future code developers to improve upon the analytical work presented herein and by 
designers to help in the preliminary sizing of friction dampers for turbine blade 
applications. 
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ABSTRACT 


The Airfoil Vibration Damper program has consisted of an analysis phase and a testing 
phase. During the analysis phase, a state-of-the-art computer code was developed, which 
can be used to guide designers in the placement and sizing of friction dampers. The use of 
this computer code was demonstrated by performing representative analyses on turbine 
blades from the Hig h Pressure Oxidizer Turbopump (HPOTP) and High Pressure Fuel 
Turbopump (HPFTP) of the Space Shuttle Main Engine (SSME). The testing phase of the 
program consisted of performing friction damping tests on two different cantilever beams. 
Data from these tests provided an empirical check on the accuracy of the computer code 
developed in the analysis phase. 

Results of the analysis and testing showed that the computer code can accurately predict the 
performance of friction dampers. In addition, a valuable set of friction damping data was 
generated, which can be used to aid in the design of friction dampers, as well as provide 
benchmark test cases for future code developers. 
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1.0 INTRODUCTION AND SUMMARY 


Turbine blades are subjected to severe loading conditions during normal operation, 
consisting of a combination of thermal, centrifugal, power bending, and oscillatory forces. 
The oscillatory forces are caused by disturbances in the hot-gas flow due to upstream and 
downstream obstructions, which generally occur at some multiple of pump speed. Dynamic 
excitation of the turbine blade occurs, as it passes through these flow disturbances during 
rotation of the disk. Although it is desirable to design turbine blades to avoid coincidence of 
natural frequencies and flow stream excitation frequencies, aerodynamic constraints 
imposed on the blade airfoil make it nearly impossible to design a completely resonant-free 
blade. Natural frequencies exist in almost all blade applications on the SSME. These 
natural frequencies are in close proximity to excitation frequencies in the hot-gas stream. 
Therefore, the turbine blade designer is forced to rely on additional damping to improve 
turbine blade fatigue life. Usually this takes the form of friction dampers, consisting of 
small, centrifugally-loaded, metal plates. These plates connect each blade to the two 
adjacent blades and provide damping, as vibration occurs, by frictional scrubbing. 

Historically, friction dampers have been designed using empirical data obtained from spin 
rig and engine testing. Analytical methods of determining friction damper performance have 
been developed, but these have proven to be either too expensive (transient time-history 
methods) or too simplistic (single-mode, steady- state harmonic balance methods) for 
practical use. Therefore, the computer code BLDAMP was developed to economically 
predict turbine blade friction damper performance. The code is not limited by the size of the 
finite-element representation of the blade. It includes the effects of as many as eight modes 
of vibration, as well as allowing input of any arbitrary forcing function. 

In addition to the analytical task described above, laboratory testing was completed to 
determine the effectiveness of friction damping in the reduction of vibratory stresses in a 
cantilever beam. This testing was done to provide a good empirical database with which to 
validate the computer code and also to serve as a benchmark for future code developers and 
damper designers. The tests were performed using two different sized beams; one with a 
low natural frequency and one with a high natural frequency. Two different dampers of the 
same geometry, but different materials, were tested on each of the beams to determine their 
influence on beam response. Results showed that friction damping can be extremely 
effective in reducing the response levels of vibrating structures. Friction damping was 
found to provide the greatest reduction in dynamic response when the damper was located 
near the tip of the beam and also when the normal load was above a critical value. 

The friction damping test results were compared with analytical predictions made using the 
computer program BLDAMP, which was developed in an earlier phase of this program. 
The comparison showed good agreement between test and analysis results. However, the 
analytical predictions were highly dependent on friction coefficient. Since friction 
coefficients were not readily available for the materials used in the test program, this 
parameter was used to tune the analytical results to match test data. 
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2.0 TECHNICAL DISCUSSION 


2.1 BLDAMP COMPUTER CODE 

2.1.1 Summary 

The computer code BLDAMP, developed by Griffin Consulting for the task of predicting 
the performance of friction dampers, was installed on the Rocketdyne CDC com pute r 
system. It was thoroughly exercised, using the High-pressure Fuel Turbopump (HPFTP) 
first-stage turbine blade of the Space Shuttle Main Engine (SSME) as a sample case. 
Several modifications were made to improve ease of use and output format, as part of the 
code implementation. The new user manual, as well as a complete listing of the modified 
code, is included in this report. 

Computer predictions for the performance of the HPFTP first-stage turbine blade dampers 
are included as a sample exercise to demonstrate the various options of BLDAMP. The 
results presented herein should not be used for design purposes, as the forcing functions 
used in the analysis were rough assumptions, which were not based on engine 
measurements or hydrodynamic analysis of the flow field surrounding the blades. 

2.1 .2 Code Description and Use. BLDAMP is a special purpose computer code 
developed by Griffin Consulting. It analytically predicts the performance of turbine blade 
friction dampers. The theoretical development of the method used by the code is reproduced 
from portions of Reference 1 and is included as an appendix (Section 6. 1) of this report. 

Input to the code consists of modal parameters, which describe the dynamic characteristics 
of a turbine blade. These parameters are derived either from an analytical representation of 
the blade or from a modal test. Damper stiffness properties and friction data are also input 
as problem parameters. A complete description of all the required input data for running 
BLDAMP, as well as a listing of the code, is provided in Section 6.2 of the Appendix. 

The code is capable of performing three different types of computer analyses. Each analysis 
type is described briefly below. Sample cases are also presented to illustrate the option. 
The sample cases were provided by the code developer and were subsequently rerun with 
identical results, using the Rocketdyne version of the code. 

2.1. 2.1 Frequency Response Analysis. This option computes the peak 
response of a friction damper, due to a sinusoidal forcing function over a frequency range of 
interest for a given normal load. The response is printed in a table of amplitude (stress or 
displacement) versus frequency. A sample plot of the output data from the program is 
presented in Figure 2.1-1 for a two degree of freedom oscillator. It is important to 
understand that die damper normal load, for which the curve of Figure 2.1-1 is calculated, is 
actually half of the centrifugal load for a single damper. This distinction is made because, 
for each damper, half of the load is supported by each blade. The program still considers 
two dampers attached to the blade. However, both the input and output data refer to a 
normal force which is only half of the centrifugal load of one damper. 

2.1 .2.2 O ptimum Normal Force Analysis. The optimum normal force 
analysis option computes the resonant response of a blade under sinusoidal loading at 
various damper normal force values. Printed output consists of a table of peak response 
versus normal force. Thus, an optimum damper normal force can be obtained using this 
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Figure 2 . 1-1 Two degree of freedom frequency response plot 
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option. An optimum normal force plot is shown in Figure 2.1-2 for the two degree of 
freedom oscillator. Again, note that the normal force values given in the curve are equal to 
one-half of the centrifugal load due to a single damper. 

2.1 .2.3 Damper Performance Analysis. The damper performance 
analysis option determines the response of a frictionally-damped turbine blade, as a function 
of the response of the same blade with no friction damper. It is easy to determine the 
improvement gained by the addition of friction dampers from information of this type. A 
damper performance curve for a two degree of freedom oscillator is presented in Figure 
2.1-3. It clearly shows the effectiveness of friction damping. 

2.1.3 HPFTP Damper Optimization. An SSME blade was used as a test case to 
exercise the code and to gain experience with systems having the same dynamic 
ch aracte ristics as actual rocket engine turbine blades. The blade chosen for analysis was the 
HPFTP first- stage turbine blade (P/N R00 19821). Modal results from an existing finite- 
element model of the blade were available for use. The only unavailable information was 
the modal stress distribution, which was not calculated during the initial analysis. In lieu of 
rerunning the original analysis to obtain modal stresses, it was determined that blade tip 
displacement would be an adequate measure of blade response to be used in the calculations. 
Normally, an equivalent alternating stress at some critical location is used as a tracking 
parameter and indicator of the level of response. At a later date, modal stress information 
became available. The displacement-based optimization curves were used, in conjunction 
with this new i nformation, to obtain rough estimates of the optimum damper normal force 
for the HPFTP turbine blade. 

Geometry plots of the blade model used in the analysis are shown in Figure 2.1-4. 
Boundary conditions consisted of restraining the model along two rows of nodes at the 
shank root, while material properties used in the analysis reflect those at SSME operational 
temperature. Comparison with blade-in-block holographic modal test results shows an 
average error of only 6.3 percent for the first 6 modes. Centrifugal stiffening effects were 
not included in the analysis. Modal information was derived from the model and was used 
as input to BLDAMP. A Campbell diagram based on the computer-generated results from 
the model is presented in Figure 2.1-5. This Campbell di agram should not be used for 
design purposes. The official Campbell diagram for the HPFTP first-stage blade is based 
on spin test results. It correctly includes the effects of centrifugal stiffening and disk/fir tree 
interface flexibility. The analysis presented herein does not attempt to address these issues, 
but is intended to serve only as an illustrative example on the use and capabilities of the 
computer program BLDAMP. 

The frequency response option of BLDAMP was exercised by forcing the blade second 
mode with a 13N sinusoidal forcing function. Frequencies ranged from approximately 
6,000 to 12,000 Hz. The forcing function was chosen such that the generalized force for 
the second mode was of such a magnitude as to cause a unit displacement of the blade tip on 
an undamped blade being forced at resonance. All other generalized forces were chosen to 
be zero. Note that undamped refers to no friction damping, in the context of this report A 
viscous damping ratio of 0.01 was used in all cases, both with and without friction damping 
present. Figure 2.1-6 presents a comparison of the frequency response plots for the 
undamped- and frictionally-damped blades. The resonant response of the blade second 
mode, which occurs at approximately 10,100 Hz, is clearly shown in the figure. The 
frequency shifts to a higher value with the addition of friction damping to the system. This 
expected result is opposite to the effect of viscous damping, which reduces the resonant 
frequency as damping is increased. 
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Figure 2.1-4 HPFTP first stage turbine blade finite element model 
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Figure 2.1-5 HPFTP first stage turbine blade Campbell diagram 
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Figure 2.1-6 HPFTP first stage turbine blade freguency response 




The optimum normal force option of the code was exercised to a greater extent than any 
other option, during the study of the HPFTP blade. Four different resonant points were 
investigated. They are shown by the circled areas on the Campbell diagram of Figure 2.1-5. 
Recall from previous discussion that these points do not represent actual engine interference 
points, because the blade model reflects the root-fixed condition. During this phase of the 
study, the damper stiffness and coefficient of friction were varied to determine die effect on 
optimum normal load. Results are presented in Figures 2.1-7 through 2.1-10. They reveal 
the classic damper curve, clearly showing the desired optimum damper normal force. These 
figures are normalized, such that the response of the blade tip is 1.0 inches under a no- 
damper configuration. Scaling to another displacement and normal load can be done in a 
linear fashion. For example, consider the second blade mode being excited by 13N. 
Consulting Figure 2.1-7 at a damper stiffness of 500,000 lb/in and an undamped tip 
displacement of 0.001 inches, the optimum 1/2 normal load for this mode is approximately 
13 pounds, whic h gives an optimum total damper load of 26 pounds. For comparison, the 
current HPFTP damper normal force is shown in Figure 2.1-11, as a function of pump 
speed. It is 210 pounds at 36,000 rpm. However, this does not imply that the current 
damper is too heavy, because of the initial assumption on the magnitude of the forcing 
function. If an undamped tip deflection twice as large had been assumed, then the optimum 
normal force would also be doubled. It is also apparent from Figures 2.1-7 through 2.1-10 
that the optimum normal force varies, depending on the mode which is being damped. 

Another set of curves can be generated, from the families shown in Figures 2.1-7 through 
2.1-10, by considering only the optimum points and by choosing damper stiffness as the 
independent variable. These are shown in Figures 2.1-12 through 2.1-15. Here, the 
generalized force is scaled such that it drives the undamped blade to a peak tip deflection of 
0.001 inches at resonance. These curves are highly dependent on the value of force used. 
For example, if the driving force magnitude is great enough to cause a tip displacement of 
0.002 inches in an undamped blade, then the optimum damper loads shown in the figures 
will double. Therefore, it is clear that, when using the optimum damper force approach, the 
forcing function magnitude must be known to a relatively high degree of accuracy to obtain 
meaningful results. 

As mentioned previously, modal stresses for the blade became available after the curves of 
Figures 2.1-7 through 2.1-15 were generated. The maximum alternating stress was then 
determined from the modal stress information for each of the first six modes, assuming a tip 
deflection of 0.001 inches. The results for each of the four interference points were then 
scaled to an undamped stress level of 15,000 psi, which was assumed to be representative 
of an undamped turbine blade stress level. The results of this scaling are presented in 
Figures 2.1-16 through 2.1-19, which represent the best estimate of the optimum damper 
weight for the HPFTP first-stage blade. 

The damper performance curve option of the code was used on a limited basis during this 
study. Results of several cases, ran for the second mode under 13N excitation, are 
presented in Figure 2.1-20. They show the relative merits of the frictionally-damped blade, 
as compared to the undamped blade. 

2.1.4 Com ments on the Use of BLDAMP. The computer code BLDAMP 
represents a significant improvement in capability over previous methods of analyzing 
friction dampers. An estimate of the optimum friction damper size, for a given blade 
design, can be made in a very short period of time, using this code. However, the code is 
not without its problems. BLDAMP is very sensitive to the input data for the analysis. 
Small input changes can make the difference in obtaining a correct or incorrect answer. 
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Figure 2.1-8 HPFTP turbine blade damper optimization curve for mode 3 
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Figure 2.1-9 HPFTP turbine blade damper optimization curve for mode 4 
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Figure 2.1-11 HPFTP turbine blade damper normal force 
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In particular BLDAMP fails to converge to a soluti on at higher values of damper stiffness. 
As an example, the damper stiffness for the HPFTP blade has been calculated to be 
1,000,000 lb/in. Yet, the program will not converge at stiffnesses larger than 200,000 lb/in 
for the fourth mode (See Figures 2.1-9 and 2.1-14). In other cases, the program will 
converge to an answer that is obviously incorrect. This can be observed from Figure 2.1- 
10, where the k = 500,000 lb/in case does not give reasonable results. Although BLDAMP 
is a great improvement over previous methods of damper analysis, additional work is 
needed to improve the convergence properties of the program. 

Other problems, which were discovered when using BLDAMP during this study, are listed 
below. 


1. The analysis does not converge at intermediate values of damper stiffness and 
then does converge at higher stiffnesses. 

2. Sometimes the program prints too few points to define the normal force curve. 

It needs user control of the normal force spacing. 

3. The analysis does not converge for low values of input force. Problems are 
anticipated when using realistic forcing functions. 

2.1.5 Conclusions. The computer code BLDAMP has been modified and installed on 
the Rocketdyne CDC computer system. Sample cases have been run, which exactly match 
those provided by the code d evelo per, Griffin Consulting. In addition, more realistic cases 
have been run, using the HPFTP first-stage turbine blade as a test case. Results of the 
study show that the program is easy to use and provides answers which appear to be 
reasonable. Based on an assumed forcing function, which gives an alternating stress of 
15,000 psi on an undamped blade, the optimum damper weight was calculated to be 
approximately 140, 10, 25, and 100 lbs for modes 2, 3, 4, and 6, respectively. All of these 
values are less than the current 210 lb damper weight. However, these damper weights are 
to be used only for illustrative purposes, because the forcing function used for the analysis 
is not accurate. 
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2.2 FRICTION DAMPING TEST RESULTS 
2.2.1 Initial Testing 


2. 2. 1.1 Summary. A laboratory test was completed to determine the 
effectiveness of friction damping in the reduction of vibratory stresses in a cantilever beam. 
The test was performed using two different sized beams; one with a low natural frequency 
and one with a high natural frequency. Two different dampers of the same geometry, but 
different materials, were tested on each of the beams to determine their influence on beam 
response. Results showed that friction damping can be extremely effective in reducing the 
response levels of vibrating structures. Friction damping was found to provide the greatest 
reduction in dynamic response when the damper was located near the tip of the beam and 
also when the normal load was above a critical value. 

2.2.1 .2 Test Objectives. The objectives of the vibration test were to measure 
the performance of friction dampers and to quantify certain parameters that were found to be 
important. These parameters are friction coefficient, damper stiffness, damper location, and 
beam driving force. The test was also intended to serve as a database for comparison with 
analytical results, obtained using the BLDAMP computer code. 

2.2.1 .3 Test Hardware and Setup. The test was repeated on two entirely 
different cantilever beams; one which was designed to operate at relatively low frequencies, 
ranging from 300 Hz to 500 Hz, and another which was designed for the 2,000-3,000 Hz 
frequency range. The low-frequency test beam (long beam) shown in Figure 2.2-1, had a 
measured first natural frequency of 344 Hz and a second natural frequency of 1,790 Hz, 
with no dampers installed. The large natural frequency spacing ensured no coupling 
between modes. The high-frequency test beam (short beam), shown in Figure 2.2-2, had 
measured undamped natural frequencies of 2,556 Hz and 4,082 Hz. 

During testing, the friction damper was installed in a selected position underneath one of the 
tangs, which were machined as integral parts of each beam. Figure 2.2-3 shows a damper 
as installed in one of the test beams and Figure 2.2-4 shows the dimensions of the damper, 
of which several samples were constructed of both Haynes 188 and silicon nitride. Only 
one of the two material types and only one damper were used on any single test run. The 
choice of Haynes 188 was made because this has been used successfully at Rocketdyne as a 
damper material on both the high-pressure oxidizer and high-pressure fuel turbopumps. 
Silicon nitride was chosen as an experimental material for its low coefficient of friction, low 
density, and high stiffness. These properties are all desirable in friction dampers, because 
past experience has shown that heavy dampers, or dampers with high friction coefficients, 
are subject to lockup. Also damper stiffness has been shown analytically to provide better 
damper performance. 

The dampers were designed to rest in a groove machined into the stationary part of the test 
fixture (Figures 2.2-1 through 2.2-3). The knife-edge portion of the damper was designed 
to fit tightly into the groove, thus precluding any damper motion relative to the base. This 
design allows only sliding motion at the interface between the damper and the test beam. 
The normal force was applied to the back side of the damper by means of a cable or string 
arrangement, which acted through a damper holder. The damper holder, shown in Figure 
2.2-5, was identical for low- and high-frequency beams. It served to transfer the load from 
the cable to the damper in a uniform fashion. 
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Figure 2.2-1 Low frequency test beam, support, and spacer 
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Figure 2.2-2 High frequency test beam, support, and spacer 
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Figure 2.2-3 Details of damper installation and loading method 
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Figure 2.2-4 Damper for high and low frequency beams 
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Figure 2.2-5 Damper holder for hicrh and low frequency beams 
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The test setup consisted of a cantilevered beam bolted to a shaker table, as shown in Figures 
2.2-6 and 2.2-7. The low-frequency beam was tested using a MB C-10 shaker in EDL 
(Engineering Development Laboratory) Vibration Room 3. Later, a Ling A249 shaker in 
Vibration Room 4 was used. The high-frequency beam was tested using a Wilcox Research 
D125L-M shaker. As can be seen from the figures, the test setup had provisions for placing 
a friction damper at any one of several locations along the span of the beam. Dampers were 
loaded by means of a cable and load cell, or by a string and pulley arrangement, as shown in 
Figures 2.2-8 and 2.2-9. The load cell was used primarily for the higher loads, while the 
string and pulley were used for the lower loads. Loads were measured either by the load 
cell or by a small spring scale. Application of the loads was accomplished by addition of 
weights to a small bucket, or by tightening a bolt attached to the load cell. 

2^2.1 .4 Test Ins trumentation. Each beam was instrumented with strain 
gages and accelerometers, positioned as shown in Figures 2.2-10 and 2.2-11. All data 
channels were recorded on magnetic tape during the test . A schematic of the instrumentation 
and data-gathering system is shown in Figure 2.2-12. 

2*2/1 -5 Modal Testing. Both high- and low-frequency beams were modal 
tested to determine natural frequencies and mode shapes. Testing was accomplished with 
the HP Modal Analyzer system, using the calibrated hammer method. Measurements were 
taken with and without a damper installed, to determine the range of resonant frequencies to 
expect when performing the actual damper sine-sweep testing. The data were also used to 
determine the damper stiffness, as will be discussed later. When the damper was installed, 
a large normal force was used to keep the damper from slipping during the test. This was 
done to ensure that the test with the damper installed represented an upper bound to the 
frequency which could occur. Table 2.2-1 presents a list of the modal frequencies found 
from the test. Figures 2.2-13 through 2.2-17 show typical mode shapes and frequency 
response functions of the beams, both with and without the damper installed. All mode 
shapes, measured in the test, consisted of beam bending in the vertical plane. 

2.2.1. 6 Dampe r Stiffness. Damper stiffness was determined by measuring 
the natural frequency of the beam, both with and without the damper installed Then a finite- 
element model of the beam was used to determine the additional stiffness necessary to cause 
the observed frequency shift. The finite-element model was calibrated by comparing 
analytical frequencies with measured frequencies, without the damper installed. The model 
was modified until a satisfactory match was obtained. When the beam was tested with the 
damper installed, a large normal force was used to ensure that the damper was completely 
locked. This had the effect of placing a spring to ground at the damper position. The 
stiffness of the damper spring was found by using the finite-element model to recompute the 
natural frequencies of the beam, as a function of damper stiffness. Then the stiffness was 
determined which gave the correct frequency, as measured in the test. This procedure was 
repeated at each damper location. An average stiffness of 450,000 Ib/in resulted for the 
Haynes 188 damper. The stiffness of the silicon nitride damper was determined from the 
Haynes damper stiffness, by using the modulus ratio between the two materials. This 
resulted in a stiffness for the silicon nitride damper of 675,000 Ib/in. 

. 2.2.1 .7 Da mper Effectiveness Testing. Testing to determine the 
effectiveness of friction dampers was accomplished by performing a sine sweep through the 
frequency range of the first beam bending mode and recording the accelerometer and strain 
gage responses. The tests were first performed without the damper installed and then with 
the damper installed. The damper normal load was selected prior to the test. Then it was 
held constant throughout the sine sweep. For the next sweep, a new value of damper 
normal force was used. The test was then repeated at the same input level. Typical 
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Table 2.2-1 

BEAM NATURAL FREQUENCIES WITH AND WITHOUT DAMPERS INSTALLED 

MODAL TEST RESULTS 


LOW FREQUENCY BEAM 


CONFIGURATION 




MODE 1 

MODE 

2 

No Damper 





344 Hz 

1790 

HZ 

Haynes 

188 

Damper 

at 

Tang 

1 

405 HZ 

>2000 

HZ 

Haynes 

188 

Damper 

at 

Tang 

3 

712 HZ 

>2000 

HZ 

Haynes 

188 

Damper 

at 

Tang 

5 

1288 Hz 

1605 

Hz 

Haynes 

188 

Damper 

at 

Tang 

7 

1416 Hz 

>2000 

HZ 

HIGH FREQUENCY 

BEAM 







CONFIGURATION 




MODE 1 

MODE 

2 

No Damper 





2556 HZ 

4082 

HZ 

Haynes 

188 

Damper 

at 

Tang 

1 

2633 Hz 

4402 

Hz 

Haynes 

188 

Damper 

at 

Tang 

2 

4335 HZ 

6052 

HZ 


Notes: (1) Tang 5 location was very close to the nodal point for the 
low frequency beam 2nd bending mode. 

(2) When dampers were installed a large normal force was applied 
to lock the damper and prevent sliding during the modal 
test. Testing was done using very light hammer raps to 
excite the beams. 

(3) All modes are beam bending modes. 
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Figure 2.2-9 High frequency beam damper loading mechanism 
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Figure 2.2-10 Low frequency beam strain gage instrumentation 
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Figure 2.2-11 High frequency beam strain gage instrumentation 
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Figure 2.2-12 Data recording system 
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Figure 2.2-13 Frequency response function and mode shape of 
low frequency beam without damper 
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Figure 2.2-14 Frequency response function and mode shape of 

low frequency beam with locked damper at tang 1 
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Figure 2.2-15 Frequency response function and mode shape of 
low frequency beam with locked damper at tang 
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response plots are shown in Figures 2.2-18 through 2.2-27 for both damped and undamped 
beams. The input was increased after testing was completed at a given excitation level, so 
that the effects at higher vibration levels could be determined. 

2.2.1 .8 Test Results. Results of the testing were plotted in the form of stress 
versus normal force curves for each damper position and type. These curves are presented 
in Figures 2.2-28 through 2.2-37 for the low-frequency beam (long beam) and Figures 2.2- 
38 through 2.2-44 for the high-frequency beam (short beam). The stress was determined 
from strain gages located at die root of the beam for all the cases shown. The normal force, 
presented in these curves, represents the total force applied to the damper by the cable. The 
actual force, at the interface between the damper and vibrating beam, is half of the total 
force. The curves show the classic friction damper behavior, which consists of a reduction 
in vibratory response as normal force is increased up to a critical value. References 5 and 6 
show a large increase in response at normal loads greater than the critical value. However, 
in this test, a flat or slight increase in response was observed as normal force was increased. 
This is due to the relatively high stiffness of the damper, which has the effect of flattening 
the response curve. 

Viscous damping ratios were computed for each beam, without dampers installed, using the 
half power point method. Damping ratios were also computed at each of the two vibration 
input levels used for the beams. These are shown in Table 2.2-2 and range from 0.005 to 
0.006 for the low-frequency beam and from 0.008 to 0.015 for the high-frequency beam. 

The major observation from the test data is that friction damping can be extremely effective 
in reducing vibratory response. Tenfold reductions in response were commonly observed 
during testing, as can be seen from the summary given in Table 2.2-3. In addition, it was 
found that the effectiveness of the dampers increases, as the damper is moved to locations 
closer to the tip of the beam. This expected result occurs because the beam motion is larger 
at the tip and can be expected to dissipate more energy at that location. This can be seen 
from a comparison of Figures 2.2-29 and 2.2-33, which show the response of the low- 
frequency beam with dampers at tangs 1 and 3, respectively. Also, from these two figures, 
it can be seen that the optimum normal force is significantly reduced, as the damper is 
moved from tang 1 to tang 3, i.e., closer to the tip of the beam. Further testing with the 
damper located outboard of tang 3 was attempted with mixed results. For these 
configurations, the beam response was consistently low level. However, the data quality 
and test-to-test repeatability were poor. It was found that, because of feedback caused by 
the stick-slip motion of the damper, the shaker control system was unable to keep the base 
moving in a sinusoidal fashion at the desired input level. 

The character of the damper performance curve for the silicon nitride damper was found to 
be somewhat different than for the Haynes 188 dampers, as can be seen by comparing 
Figures 2.2-28 and 2.2-29 with Figures 2.2-30 and 2.2-31. From Figures 2.2-28 and 2.2- 
30 it is clear that the silicon nitride damper requires a much larger normal load to reach the 
optimum point than the Haynes 188 damper. This is indicative of a low friction coefficient, 
which is expected for this material. On the other hand, Figures 2.2-29 and 2.2-31 show the 
optimum normal force to be about equal for the two dampers. A closer evaluation of the 
curve of Figure 2.2-31 leads one to believe that the last two data points are in error and that 
the actual optimum point is much farther to the right than shown. This is supported by the 
fact that at 2G input (Figure 2.2-30) the optimum normal force was found to be 250 lbs. 
Therefore, at 5G input (Figure 2.2-31) the optimum must be greater than 250 lbs. 

Testing on the low-frequency beam at the outboard damper locations and on the high- 
frequency beam at all damper locations proved to be very difficult, because of damper 
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Table 2.2-2 


CRITICAL DAMPING RATIOS WITHOUT FRICTION DAMPERS INSTALLED 


TOW FREQUENCY BEAM 

EXCITATION LEVEL 
2G 
5G 


DAMPING RATIO 
. 005 
.006 


HIGH FREQUENCY BEAM 

EXCITATION LEVEL 
8G 
12G 


DAMPING RATIO 
.015 
.008 
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Table 2.2-3 

PEAK RESPONSES WITH AND WITHOUT DAMPERS 



INPUT 

STRESS 

STRESS 


(N=0) 

(N=Optimum) 

2G 

11.7 

ksi 

1.2 

ksi 

5G 

23.2 

ksi 

3.8 

ksi 

2G 

11.9 

ksi 

1.4 

ksi 

5G 

23.2 

ksi 

7.4 

ksi 

2G 

11.7 

ksi 

0.2 

ksi 

5G 

23.2 

ksi 

0.6 

ksi 

2G 

11.9 

ksi 

1.4 

ksi 

5G 

23.2 

ksi 

— 


2G 

11.7 

ksi 

— 


2G 

11.7 

ksi 

— 



N 

TANG 

DAMPER 

70 lbs 

1 

Haynes 

120 lbs 

1 

Haynes 

250 lbs 

1 

Si/Nitride 

1200 lbs 

1 

Si/Nitride 

12 lbs 

3 

Haynes 

25 lbs 

3 

Haynes 

18 lbs 

3 

Si/Nitride 

— 

3 

Si/Nitride 

— 

5 

Haynes 

— 

5 

Si/Nitride 


HIGH FREQUENCY BEAM 


INPUT 

STRESS 

STRESS 


(N=0) 

(N=Optimim) 

8G 

2.8 ksi 

0.2 ksi 

12G 

6.8 ksi 

0.8 ksi 

8G 

2.9 ksi 

0.9 ksi 

12G 

6.8 ksi 

0.7 ksi 

8G 

2 . 9 ksi 

0.5 ksi 

12G 

6.1 ksi 

0.6 ksi 

12G 

6.9 ksi 

0.6 ksi 


N TANG DAMPER 


60 

lbs 

1 

Haynes 

50 

lbs 

1 

Haynes 

75 

lbs 

1 

Si/Nitride 

50 

lbs 

1 

Si/Nitride 

50 

lbs 

2 

Si/Nitride 

45 

lbs 

2 

Haynes 

35 

lbs 

2 

Si/Nitride 
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LONG BEAM DAMPING TEST/ANALYSIS 


2G INPUT— SILICON DAMPER AT TANG1 



NORMAL FORCE (LBS) 

□ TEST DATA + BLDAMP PROGRAM yU=Q5 

Figure 2.2-30 Low frequency beam optimum normal force curve 

for 2G input and Silicon Nitride damper at tang 1 
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LONG BEAM DAMPING TEST/ANALYSIS 



0 2 4 6 8 10 12 14 16 18 20 

NORMAL FORCE (LBS) 

□ TEST DATA + BLDAMP PROGRAM ^J=.2 


Figure 2.2-32 Low frequency beam optimum normal force curve 
for 2G input and Haynes 188 damper at tang 3 
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Figure 2.2-34 Low frequency beam optimum normal force curve 

for 2G input and Silicon Nitride damper at tang 3 
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Figure 2.2-35 Low frequency beam optimum normal force curve 
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LONG BEAM DAMPING TEST/ANALYSIS 
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Figure 2.2-37 Low frequency beam optimum normal force curve 
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Figure 2.2-38 High frequency beam optimum normal force curve 
for 8G input and Haynes 188 damper at tang 1 


TEST DATA 



— 

m w 

CL v 

^ c 


in J5 
m a 

UJ o 

(£ .c 



to 


o 
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SHORT BEAM DAMPING TEST/ANALYSIS 
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□ TEST DATA + BLDAMP PROGRAM ^=05 

Figure 2.2-40 High frequency beam optimum normal force curve 

for 8G input and Silicon Nitride damper at tang 1 
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SHORT BEAM DAMPING TEST/ANALYSIS 


1 2G INPUT - SILICON DAMPER AT TANG 1 
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Fiaure 2.2-41 Hiqh frequency beam optimum normal force curve 
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Figure 2.2-42 High frequency beam optimum normal force curve 

for 8 G input and Haynes 188 damper at tang 2 
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chatter and shaker control problems. As the damper transitioned from a locked state to a 
slipping state, a loud noise emanated from the damper area and the shaker could no longer 
be controlled. The damper would not remain centered in position and tended to drift off to 
one side. It was also found that the output signals from the strain gages and accelerometers 
were not sinusoidal during this period. The nonsinusoidal nature of the vibration is 
significant, because the BLDAMP program assumes sinusoidal excitation. 

The surface finish on the test beams was originally machined smooth. However, after a 
short period of testing, rub marks were observed at the damper contact points. These rub 
marks were very similar to the marks observed on the HPOTP blade platforms after time in 
service. The contact surface for the Haynes 188 dampers, however, was worn in a much 
different manner than dampers that have been hot-fire tested on actual blades. The damper 
surface was extremely rough and pitted. A black powdery residue was observed on both 
the beam and damper after testing. The residue is believed to be coming from the damper, 
because most of the wear occurred there. The reason behind the abnormal wear patterns is 
most likely because of the dynamic response of the damper itself. As will be discussed later 
in this report, die damper was vibrating excessively during the period of time when the 
beam was passing through resonance. The silicon nitride damper contact surface did not 
show the wear that was present on the Haynes damper. 

2.2.1 .9 Problems. As can be seen from a comparison of the plots of Figures 
2.2-23 and 2.2-25, the character of the vibration changed dramatically when friction 
dampers were installed. The plots of these figures represent RMS data and do not show the 
actual time histories of the signal. Although the time histories are not presented here, the 
intent was to keep the input to the shaker as sinusoidal as possible, during the entire sweep 
through the frequency range of interest. This proved quite difficult when friction dampers 
were installed, particularly in the high-frequency beam and also in the low-frequency beam 
at damper locations near the tip. During testing under these conditions, the beam response 
was not sinusoidal and feedback into the shaker table distorted the input, to a large degree. 
Since the analytical program BLDAMP assumes a constant magnitude sinusoidal input, it 
was impossible to correlate the analysis results to the tests, during which damper chatter 
occurred. 


2.2.2 Follow-on Testing 

2.2.2.1 Summary. Another test series was completed using the low-frequency 
beam, due to the problems encountered with the initial testing, as described in Section 
2.2.1. Testing was performed at lower input levels, using a different shaker and control 
system, which helped prevent the nonsinusoidal oscillations observed with the previous 
setup. This testing did not alter the basic conclusions regarding friction dampers; but the 
quality of the test results was improved significantly. A large database was also acquired, 
which will aid in the development of future analytical methods of damper design. 
Furthermore, a set of design curves was generated, which will greatly aid in the selection of 
the correct damper design parameters for a given blade configuration. 

2.2.2.2 Test Description. Testing was performed in an identical manner to 
that described in Section 2.2.1, with several exceptions. Most importantly, the input 
vibration levels were reduced from the previous testing. This was done because excess 
motion of the test beam at high vibration levels caused the damper to chatter. This behavior 
is probably not indicative of actual turbine damper behavior, because turbine damper blade 
motion is an order of magnitude smaller than what was used on the beam during the initial 
testing. Also, very large centrifugal forces, on the order of 150,000 Gs, are present in 
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typical turbines. These forces, which are impossible to duplicate in a nonrotating test setup, 
should keep the damper from chattering by holding it tightly against the blade surface. The 
second improvement to the test was the use of a LTV Model 275 A vibration shaker, which 
was rated at 10,000 lb force. This shaker was coupled with an Unholtz Dickie SP-7 control 
system, which helped to reduce feedback and kept the shaker motion as sinusoidal as 
possible. 

Data was obtained by first placing the damper at one of seven locations on the test beam and 
then slowly sweeping through the resonant point. Two series of tests were performed on 
the beam for each of the seven damper locations; one at a shaker input level of 0.5G peak 
and another at a level of 1.0G peak. The normal force on the damper was changed, after 
each test run, to obtain curves of response, as a function of normal force. Testing 
proceeded with both input levels and two damper types, until data were obtained at all seven 
damper locations. Over 100 test runs were completed during the course of the testing. 

2.2.2.3 Test Results. Curves of peak beam response at the resonant point, as 
a function of damper normal load, are presented in Figures 2.2-45 through 2.2-72 for the 
Haynes 188 damper and in Figures 2.2-73 through 2.2-100 for the silicon nitride damper. 
Each of the curves shows the characteristic damper behavior. This behavior consists of an 
initial drop in beam motion, until an optimum normal load is reached; at which point the 
response flattens. Further increases in normal load fail to reduce the beam response and, in 
some cases, may even increase the response slightly. 

Friction damping was shown to be extremely effective in reducing vibratory response, as 
was observed in the previous testing on this beam. Figures 2.2-101 through 2.2-106 
summarize damper effectiveness, as a function of location on the beam. These data show 
that reductions of over 95 percent can be achieved with proper placement. These curves 
have been completely nondimensionalized, thereby providing a valuable design tool for 
determining the best position for a damper. These plots show the damper should be located 
at a blade span of at least 30 percent to be most effective. The results shown in these figures 
are also consistent with previous work; in that the damper performs better as it is moved 
toward the beam tip. Vibratory motion is larger toward the tip. Therefore, more energy 
dissipation will occur with the damper at that location. However, the curves are relatively 
flat, once the 30-40 percent span position is reached, which was an unexpected result from 
testing. This is an indication that a damper placed at 40 percent span is nearly as good as a 
tip damper. Tip dampers have been shown in spin testing to perform extremely well (see 
Reference 1). 

The normal force required to attain minimum response is presented in Figure 2.2-107 for the 
first mode of the beam. This curve indicates that a lower value of force is needed, as the 
damper is moved closer to the tip. Values plotted in this figure have also been 
nondimensionalized, which will aid future designers in the selection of the optimum normal 
force for their application. Proper use of this curve involves specification of the turbine 
blade generalized force, which, in turn, requires knowledge of the flow field forcing 
function and the blade mode shape. Blade mode shapes different from those of a cantilever 
beam will invalidate the results. Nevertheless, the curve of Figure 2.2-107 yields a starting 
point for selection of the optimum damper weight Also worth noting is the close overlay of 
the Haynes 188 and silicon nitride damper curves. This indicates the two damper materials 
have roughly the same friction coefficient. Friction coefficients of materials different from 
Haynes 188 or silicon nitride will result in different curves. However, most materials used 
for dampers will probably have about the same friction coefficient indicating that the curves 
will still be useful. 
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Figure 2 „ 2-45 Low frequency beam optimization curve based on test data 
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Figure 2.2-46 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 1, 0.5G input 
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Figure 2 . 2-47 Low frequency beam optimization curve based on test data 
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Figure 2.2-48 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 2, 0.5G input 
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Figure 2=2-49 Low frequency beam optimization curve based on test data 




0£3*I6 OH/ra 
Strain (u in/in peak) 


Friction Damper Performance 


T ang 3 Haynes 186 Damper 0.5G Input 



Normal Load 

□ Root S/G a Mid-Span S/G 


Figure 2.2-50 Low frequency beam optimization curve based on test data 

Haynes 188 damper at tang 3, 0 „ 5G input 
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Fiaure 2 „ 2-51 Low frequency beam optimization curve based on test data 
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Figure 2.2-52 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 4 e 0.5 G input 
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Fiqure 2.2-53 Low frequency beam optimization curve based on test data. 
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Figure 2.2-54 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 5, 0=5G input 
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Figure 2.2-56 Low frequency beam optimization curve based on test data 

Haynes 188 damper at tang 6, 0.5G input 
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Figure 2.2-57 Low frequency beam optimization curve based on test data 

Havnes 188 damper at tang 7, 0 . 5G input 
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Figure 2.2-58 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 7, 0.5G input 
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Figure 2,2-59 Low frequency beam optimization curve based 
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Figure 2.2-60 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 1, 1.0G input 
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Figure 2.2-62 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 2, 1.0G input 
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Figure 2.2 63 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tana 3. - ~~ -* 
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Figure 2.2-64 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 3, 1.0G input 
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Figure 2.2-65 Low frequency beam optimization curve based on test data 

Haynes 188 damper at tang 4 , 1 „ 0G input 
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Figure 2.2-67 Low frequency beam optimization curve based on test data 

Haynes 188 damper at tang 5 1 1. OG input 
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Figure 2.2-69 Low frequency beam optimization curve based on test data 
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Figure 2.2-70 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 6, 1.0G input 
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Figure 2.2-71 Low frequency beam optimization curve based on test data 
Havnes 188 damper at tana 7. 1.0G input 
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Figure 2.2-72 Low frequency beam optimization curve based on test data 
Haynes 188 damper at tang 7, 1.0G input 
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Figure 2.2-75 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 2 , 0 „ 5G input 
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Figure 2.2-77 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 3, 0 . 5G input 
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Figure 2.2-78 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 3, 0.5G input 
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Figure 2 „ 2-79 Low frequency beam optimization curve based on test data 

silicon Nitride; damoer at tana 4, 0.5G Tnnut 
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Figure 2.2-80 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 4, 0.5G input 
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Figure 2.2-81 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 5, 0 . 5G input 
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Figure 2.2-82 
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Figure 2.2-83 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 6, 0.5G input 
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Figure 2.2-84 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 6, 0.5G input 
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Figure 2.2-85 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 7 , 0.5G i: 
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Figure 2.2-87 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 1, 1.0G input 
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Figure 2.2-88 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 1, 1.0G input 
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Figure 2.2-89 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 2 , 1 . OG input 
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Figure 2.2-90 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 2, 1.0G input 
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Figure 2.2-91 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 3 , 1 . OG input 
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Figure 2.2-92 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 3, 1.0G input 
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Figure 2.2-93 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 4 , 1.0G input 
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Figure 2.2-94 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 4, 1.0G input 
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Figure 2.2-95 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tana 5.1 
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Figure 2.2-96 Low frequency beam optimization curve based on test data 
Silicon Nitride damper at tang 5, 1.0G input 
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Figure 2.2-97 Low frequency beam optimization curve based on test data 
Silicon Nitride damner at tana fi i . no i nrm+- 
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Figure 2.3-7 Peak blade jet force in spin test for 19N 

excitation, undamped 
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Figure 2.3-6 Peak blade jet force in spin test for 12N 
excitation , undamped 
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type. A coefficient of friction of 0.2 was used, based on the test/analysis comparison of 
Section 2.3-1. This coefficient is consistent with the material and surface finish of the 
blades and dampers. 

The damper for the HPOTP first-stage blade is shown in Figure 2.3-3. It is a first- 
generation-type damper that stands vertically between the blades underneath the platforms. 
The damper is forced against the platforms of neighboring blades, with a total force of 38 
lbs, and provides damping by frictional means during operation. 

The primary modes of interest are the first two modes of the blade; the tangential mode and 
the axial mode. The natural frequency of the tangential mode, as measured in the spin test, 
is approximately 4,800 Hz. The natural frequency of the axial mode is approximately 9,500 
Hz. This compares with the tangential and axial modes of the model, which are 5,089 Hz 
and 11,345 Hz, respectively. The first mode can be excited by the 12N excitation at 
approximately 24,000 rpm or by the 19N excitation at approximately 14,000 rpm. The 
second mode can be excited by the 19N excitation at approximately 30,000 rpm. Spin test 
data for the HPOTP first-stage blades without dampers is shown in Figures 2.3-4 and 2.3- 
5. The corresponding peak jet forces are shown in Figures 2.3-6 and 2.3-7. Data from 
ramp number 4 (test 3-1A-4, 12N excitation, undamped) and ramp number 1 (test 3-1A-1A, 
19N excitation, undamped case) were used for the comparison. The corresponding friction- 
damped test data are shown in Figures 2.3-8 through 2.3-1 1. 

An initial run of the program BLDAMP was made to find the generalized force scaling 
factor, due to the estimation of the peak jet forces. This was done so that analysis results 
would match test data for the undamped case. The case of 19N excitation of the axial mode 
was used as the scaling point. A generalized force scaling factor of 3.02 was calculated for 
this case. This factor was then used for the remaining cases. A parametric study was then 
performed using various damper stiffnesses, once the generalized forces were fixed. The 
case of 12N excitation of the tangential mode was used to determine this parameter. Results 
are shown in Figure 2.3-13. The graph shows that the undamped data point (damper force 
= 0.0) agrees very well with test data. Also, the parametric data shows that a damper 
stiffness of 20,000 lb/in agrees very well with the friction-damped test data. The response 
curve in Figure 2.3-13 can be compared to the test data in Figure 2.3-8. This again shows 
the friction-damped amplitudes to txs very close. The program BLDAMP predicted a shift in 
frequency of resonant response of approximately 4 percent for this case. The test data 
shows a shift of approximately 12 percent 

The case of 19N excitation of the axial mode is shown in Figure 2.3-14. The undamped test 
data matches the analysis exactly, because the generalized force was scaled from this case. 
The friction-damped case also agrees quite well with the test data. The response curve in 
Figure 2.3-15 can be compared to the test data in Figure 2.3-9. The program BLDAMP 
predicted a shift in frequency of resonant response of approximately one percent for this 
case. The test data also shows a negligible shift in frequency. 

The case of 19N excitation of the tangential mode is show in Figure 2.3-16. The graph 
shows that the undamped data point is somewhat higher than test data. This could be 
caused by the equivalent viscous damping in the blade being higher than the one percent 
used in the analysis. This data point is taken at a much lower speed. The equivalent 
viscous damping could be higher, because the blade is not fully locked in the fir-tree area. 
However, the friction-damped response agrees well with the test data. To get a more 
realistic comparison, the case could be rerun with the generalized force adjusted, so that the 
undamped case matches test data. The response curve in Figure 2.3-17 can be compared to 
the test data in figure 2.3-9. The program BLDAMP predicted a shift in frequency of 
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Figures 2.2-40, 2.2-41, and 2.2-44 for the high-frequency beam. Results show the silicon 
nitride damper required a higher value of friction force to achieve the same reduction in 
vibration amplitude as the Haynes 188 damper. This indicates that silicon nitride has a 
lower coefficient of friction. The test/analysis match for the low-frequency beam, with the 
damper at tang 1, gives a friction coefficient of p. = 0.05 (see Figure 2.2-30). This value is 
approximately four times less than the Haynes 188 friction coefficient, under the same test 
conditions. This is in disagreement with the results predicted in Section 2 . 2 . 23 , which are 
based on data of a better quality than that used here. 

2.3.2 Spin Test Comparison 

The test results from a Rocketdyne High-pressure Oxidizer Turbopump (HPOTP) first- 
stage blade (P/N RS007707) spin test (see Reference 1) were compared to predictions from 
the program BLDAMP. The comparison showed good agreement between test and 
analysis results. However, the analytical predictions were highly dependent on damper 
stiffness. Since the HPOTP first-stage blade damper stiffness was not readily available, 
this parameter was used to tune the analytical model to match test data for one condition. 
Two other conditions, utilizing the same damper stiffness as in the first condition, also 
matched data very well. In addition, this comparison shows the degree to which the 
damper has been optimized. 

The Rocketdyne spin pit accomodates a turbomachinery shaft, disk, and turbine blades. 
This assembly is spun up to operational speeds in a vacuum by an electric motor. Various 
fluid forces can be placed on the blades by means of nitrogen gas jets. In one test series, 12 
evenly spaced jets of gas were directed onto the trailing edges of the HPOTP first-stage 
turbine blades. A schematic of the jet configuration is shown in Figure 2.3-1. In another 
test series, 19 evenly spaced jets were directed onto the blades. The blade response is 
measured by strain gages placed on the blades. The strain gage wires are attached to a slip 
ring, that allows the signal to be sent from the rotating blades to the stationary test 
equipment. 

The program BLDAMP and a HPOTP first-stage blade finite-element model were used to 
analytically predict the friction-damped response. The finite-element model is a 
STARDYNE plate model (shown in Figure 2.3-2). The turbine blade model is fixed at the 
root. Modal data from the model were used as input to the BLDAMP program. Modal 
stresses were obtained from the trailing-edge root, pressure side location (element number 
386 shown in Figure 2.3-2). This location matches the strain gage location in the spin test 
data. 

Generalized forces for each mode are also needed as input to the BLDAMP program. 
Because the blade will have peak responses only when it is being forced at one of its natural 
frequencies, only the harmonic components of force that match the natural frequencies of the 
blade need to be input to BLDAMP. First, the peak jet-pulse force on a blade was 
estimated, using hydrodynamic principles. Then, because the harmonic components of a 
pulse are dependent on jet spacing, a Fourier series decomposition was performed for the 12 
jets and 19 jets, individually. The first harmonic of the 12 jets was found to be 0.09163 for 
a unit pressure pulse, while the first harmonic for 19 jets was 0. 14383. The peak jet force is 
multiplied by the harmonic factor to give the harmonic force. The generalized forces were 
obtained by first matrix-multiplying the modal vectors with a unit vector, in the direction of 
the jet force, and then multiplying that resulting vector by the harmonic force. 

Equivalent viscous damping of 1 percent was used in the analysis to represent damping due 
to fir-tree motion and losses in the blade material. This value is typical for hardware of this 
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2.3 COMPARISON OF ANALYSIS AND TEST RESULTS 
2.3.1 Nonrotatina Beam Comparison 


The friction damping test results of Section 2.2.1 were compared with analytical predictions 
made using the the computer program BLDAMP, described in Section 2.1. This 
comparison showed good agreement between test and analytical results. However, the 
analytical predictions were highly dependent on friction coefficient. Since friction 
coefficients were not readily available for the materials used in the test program, this 
parameter was used to tune the results to match test data. 

Validation of the computer code was accomplished by creating finite-element models of the 
two beams, analytically simulating their response, and comparing analytical results to test 
data. Modal parameters of the two test beams and stiffness properties of the damper were 
input to BLDAMP to obtain the damper performance curves, which are presented in Figures 
2.2-28 through 2.2-44. These figures also contain plots of the test data, as was discussed 
previously. A sensitivity study was performed, with friction coefficient as the variable 
parameter, to obtain a good analytical match with the test data. This parameter was adjusted 
until the results from BLDAMP closely matched test results, since friction coefficients are 
not well known. This match is clearly shown in Figure 2.2-28, where the test and analytical 
results from the low-frequency beam with the Haynes 188 damper are compared, for three 
different values of friction coefficient. This figure indicates the best match was found by 
choosing |i = 0.2 as the friction coefficient. The test was then repeated at a higher G level, 
since friction damping is a nonlinear phenomenon. The test/analysis comparison was then 
made, using the same friction coefficient found in the previous lower-level test. This 
comparison is presented in Figure 2.2-29. The comparison shows a good match between 
test and analysis results. A similar procedure was performed for the high-frequency beam 
with the Haynes 188 damper, the results of which are shown in Figure 2.2-38. Here, the 
friction coefficient necessary for a good test/analysis match was p. = 0.05, which is 
unrealistically low and does not agree with the value obtained in the low-frequency beam 
tests. Subsequent testing of the same configuration at a higher vibration level gave the 
results shown in Figure 2.2-39. The coefficient of friction used in this figure is the same as 
was obtained from the low-level testing and yields a poor comparison with the test data. 
The poor comparison observed is probably due to damper chatter and nonsinusoidal motion 
of the shaker, which was prevalent on this beam. Therefore, for this reason, the test data 
for the high-frequency beam should be regarded as indicative of general trends only and 
these test data should not be used for detailed comparison with analytical results. 

A possible explanation, for the abnormally low friction coefficient found for the short beam 
(Figure 2.2-38), is that the damper may not have been in contact with the beam, during the 
entire cycle of vibration. As mentioned previously, when the beam passed through the 
resonant frequency, a loud noise emanated from the damper area and the damper tended to 
wander in a direction parallel to the damper knife-edge slot. It was difficult to keep the 
damper centered in the slot. Subsequent inspection of the dampers indicated considerable 
wear at the damper/test beam interface. The wear pattern observed on the dampers was not 
indicative of tangential sliding motion, as would be expected. The pattern suggested the 
damper was beating against the beam in a direction normal to the contact surface. It is 
believed the damper was not in contact with the beam, during the entire vibration cycle. 
Since the analysis assumes damper contact at all times, a lower coefficient of friction would 
be required to match the test data, which was taken with only partial contact. 

A comparison of the silicon nitride damper test data with analytical results is shown in 
Figures 2.2-30, 2.2-31, 2.2-34, 2.2-35, and 2.2-37 for the low-frequency beam and in 
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Figure 2=2-107 Optimum normal force as a function of damper position 
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Figure 2.2-104 Relative response of beam tip as a function of damper 
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Figure 2.2-101 Relative response of beam tip as a function of damper 

position for Haynes 188 damper 
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Figure 2.2-99 Low frequency beam optimization curve based on test data 

Silicon Nitride damper at tang 7, 1 . OG input 
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Figure 2.3-12 Optimization curve from BLDAMP, mode #1, 12N 
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Figure 2.3-15 Response curve from BLDAMP, mode #2, 19N 

excitation 
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Figure 2.3-16 Optimization curve from BLDAMP, mode #1, 19N 
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Figure 2.3-17 Response curve from BLDAMP, mode #1, 19N 

excitation 




resonant response of approximately 8 percent for this case, while the test data show a 
negligible shift in frequency. 

The optimization curves in Figures 2.3-12, 2.3-14, and 2.3-16 indicate that the damper is, 
for its stiffness, about twice as heavy as it should be for optimal damping of the first two 
modes. Also, the figures indicate that a stiffer damper could further reduce resonant 
response and decrease the optimum normal force sensitivity. 

Overall, the test and analysis correlated very well, although there are a few key parameters 
that need to be accurately predicted before the program BLDAMP can be truly useful as a 
design tool. The damper stiffness, forcing function, and coefficient of friction are very 
important in predicting friction-damped response. Comparisons of test data with analysis, 
such as the one described here, are helpful in the determination of these key parameters for 
future design work. 
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3.0 CONCLUSIONS 


A state-of-the-art computer code has been developed for the analysis of friction dampers, as 
a result of this program . The new code represents a significant improvement over existing 
analytical tools; in that multiple modes of the structure are accounted for, without the use of 
expensive time integration techniques. This allows the accurate representation of turbine 
blades, using the mode shapes generated from high-fidelity, finite-element models. 
Previous analytical methods considered only single modes of simple one degree of freedom 
systems. 

The new computer code, BLDAMP, was modified to improve input and output formats. 
Sample cases have been run, which exactly match those provided by the code developer, 
Griffin Consulting. In addition, more realistic cases have been run using the HPFTP first- 
stage turbine blade as a test case. Results of the study show that the program is easy to use 
and provides answers which appear to be reasonable. 

Further comparison of BLDAMP results, against test data from nonrotating test beams and 
spin pit data from instrumented turbine blades, has provided a good validation of the code. 
Results from the test/analysis comparison show that the normal force that causes the lowest 
response was highly dependent on the friction coefficient between the damper and beam 
surfaces. This makes prediction of damper performance very difficult, unless friction 
characteristics are known for the materials at hand. It was found, however, that once the 
friction coefficient was known, damper performance can be reliably predicted using the 
BLDAMP computer code. 

It can be concluded, from the results of the extensive testing presented in this report, that 
friction damping is an effective way to reduce the response of vibrating structures. 
Response reductions of over 95 percent, from the undamped configuration, were observed 
in the test beam. It was also found that the most effective position to locate the dampers is at 
a beam span of 30 percent, or more. The effectiveness of the friction damper remained 
relatively constant for damper locations between 30 and 100 percent span. This is a 
significant observation, because it was previously thought tip dampers were the only way to 
obtain response reductions of the magnitude observed in the testing. 

The test data were plotted in nondimensional form and a set of design curves was generated, 
which will greatly aid in the selection of damper parameters for future applications. On an 
equal level with die development of the BLDAMP computer code, these design curves are 
one of the major accomplishments of this program. They provide a much needed base from 
which to understand damper behavior. 
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4.0 RECOMMENDATIONS FOR FURTHER WORK 


The main problem encountered, in the initial friction damper testing reported herein, was the 
lack of good signal quality and repeatability from test to test. These problems were solved, 
for the most part, in the second test series, where the vibration amplitude was reduced and a 
more sophisticated test controller and shaker were used. The testing was still not as 
repeatable as desired and improvements could be made to improve the data quality. The 
problem centers on the method of applying the excitation to the beam, which was by use of 
base motion of the entire test apparatus. Thus, the beam was excited by inertial forces, 
which were dependent on the motion of the base. This development necessitated very 
accurate control of the shaker. As an improvement, it is suggested that further testing be 
done, using an exciter that applies an easily controllable force directly to the beam, without 
using base excitation. A magnetic excitation device may be the solution, because it can be 
commanded to provide a sinusoidal force and it would not involve any feedback, due to the 
response of the beam. 

All the test data gathered in this program was from nonrotating beam specimens. Although 
the data should be applicable to blades in a rotating environment, there is one major 
difference that may alter the results. The rotating environment differs from the nonrotating 
test in the method of application of the normal force. Centrifugal accelerations in typical 
turbines can be on the order of 150,000 Gs. These accelerations provide normal forces of 
the same magnitude as were used in the nonrotating test. However, the ratio of damper 
mass to normal force is many orders of magnitude smaller in the rotating test. The small 
mass/force ratio in the rotating environment makes it almost impossible for damper chatter to 
occur, while in the nonrotating test setup chatter was observed. A logical extension to the 
testing already completed would be to perform similar tests using the Rocketdyne whirligig 
spin facility. 

Another improvement, which would increase the quality of the data, involves damper wear 
as the test progresses. In the initial testing, a large amount of damper wear was observed 
which was attributed to damper chatter. A black residue formed on the damper contact 
surfaces, which undoubtedly affected the results by changing the frictional characteristics of 
the interface. Damper chatter was not observed, to any great extent, in the second phase of 
testing. However, a small amount of the residue appeared on the damper and beam during 
testing. An attempt was made to clean the residue as it appeared, but the cleaning process 
was not continuous. It is not known if the residue changed the results. Damper wear will 
occur in an actual turbine, but any residue will be removed immediately by the hot-gas 
stream. Therefore, an improvement to the testing would incorporate an air jet directed at the 
damper to remove the residue as it is produced. 

An extension to the testing, which may enhance knowledge of how friction dampers 
function, would involve instrumentation of the damper, as well as the beam. Currently, 
very simple models are used, which idealize the damper as a massless spring. This may not 
be adequate, as the damper may have a dynamic response of its own. For example, a strain 
gage mounted on the damper could measure the stretch that occurs before slipping. It may 
even determine if slip occurs instantaneously or over a period of time. 
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The test data was obtained on a uniform cantilever beam test specimen, vibrating in the first 
mode, only. A valuable extension to this work would be to repeat the testing for beam 
vibration in higher modes. 

Although some turbine blades have mode shapes that are similar to those of cantilever 
beams, there are many other instances where the lower modes consist of localized motion of 
the airfoil, with little or no motion of the blade as a whole. A test program could be 
developed that would generate a set of design curves for complex mode shapes, which are 
commonly found in turbine blades. An example, of this type of program, would explore the 
capability of friction dampers to eliminate trailing-edge flap modes, which exist in many 
applications. 
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1. INTRODUCTION 


The objective of this report is to summarize the underlying physical 
assumptions and mathematics that provide the basis for a computer code, 
BLDAMP, that may be used to calculate the forced response of frictionally 
damped turbine blades and optimize the design of friction dampers. BLDAMP 
will be especially suited for studying SSME turbo-pump vibration in that it 
will have the following unique capabilities: 

(1) . It will simulate mode shape changes induced by the friction 

constraints, thus it can be used to analyze dampers located at the 
tips of the blades. 

(2) . It will automatically analyze the damping of torsional modes, 

i.e. blade rotations will also dissipate energy. 

(3) . Stiff-wise blade motion as well as ease-wise motion will be 

included in the formulation, thus the code can be used to 
analyze stiffwise bending as in the 2nd mode of the HPOTP 
blade. 

These capabilities are not currently available in any existing computer code. 
As a result, the proposed computer code will significantly advance the 
state-of-the-art in friction damping analyses. The reasons that BLDAMP can 
analyze these types of problems is because a more general mathematical 
formulation is used that incorporates more physical features of the 
blade-damper system. Specifically, mode shape changes can be taken into 
account because the problem is formulated in terms of receptances rather 
than a specific mode of response.' Torsional modes may be analyzed because 
friction contact is assumed to occur at four points on each blade and since 
the points do not lie on the center of rotation the frictional forces produce a 
resisting moment which damps torsional modes. Stiff-wise blade motion is 
included in the formulation by considering that the blade may move in two 
independent directions at the damper contact point. The mathematical basis 
of each of these features is presented in this report. 

An important point to note is that the analysis of non-linear vibration 
problems is difficult. The more degrees-of-freedom or complexity that is 
introduced in modelling a problem the more computer time is required to 
solve it, and, in fact since you ultimately must solve systems of non-linear 
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algebraic equations, the algorithm may have difficulty converging if the 
system of equations is too large. As a result, developing a useful code 
involves making trade-offs between either having a tractable problem or 
including more complex physical features of the system. For this reason, it 
is important to initially model the problem as simply as possible and develop 
a code that works in an efficient manner. Once the baseline code is 
developed and tested then additional features can be incorporated if they are 
needed. 

This report first summarizes the mathematical formulation of the problem, 
gives the input requirements of the code and a flow chart, discusses 
potential enhancements to the code that may be incorporated at a later time, 
and then provides some concluding remarks. 
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2. FORMULATION 


2.1 RECEPTANCE FORMULATION OF BLA DE RESPONSE 

Initially, we consider the response on an individual blade. The aerodynamic 
forces that act on the blade are assumed to be harmonic (refer to Fig. 1), i.e. 

= f 0 TO e ,m W 

where bold letters indicate vectors or matrices. The friction forces are 
given by 

(f(t)= f.”' O) 

In both expressions the phase of the quantities are included by considering F 
and fj complex quanitities (physically the solution corresponds to the real 

parts of each of the quantities of interest). 

The response of the blade is given by the displacement vector 

Uj(t) » Uj e i<Bl (3) 


and by the stress vector 


c 




ico t 


(4) 


where i takes on integer values that correspond to points of interest on the 
blade. In particular, i equals 1 through 4 corresponds to the damper contact 
points as indicated in Fig. 1 and i = 5 is a typical reference point at which we 
wish to know the vibratory stress. 

The receptances are defined in terms of the response at a point of interest 
due to a sinusoidal excitation at a given point. In particular, the 
displacement at the ith point due to the external excitation is 


t 
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and the displacement at the ith point due to an unit excitation at the jth 
point is 



( 6 ) 


Then since the blade is linear, from (3) and (5) the displacement 
coefficients are given by 


w? + £vj 

j=i 


and analogously the stress coefficients by 

4 




(7) 


( 8 ) 


2.2 CALCULATING RECEPTANCES 

Many engineers are more familar with real receptances. In this approach if 
the excitation on a simple spring/mass oscillator is given by 

f(t) = (1) cos at (9) 


then the displacement would be given by 

u(t) =r c cos cot + r t sin cot (10) 


where r c and r s are referred to as the real receptances. Real receptances are 

readily related to complex receptances. For example, alternatively we could 
write (9) as 

f(t) = Real(e itDl ) (11) 
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thus the complex receptance, r c -ir s , is just a linear combination of the real 
receptances. 

Receptances can be calculated from modal information in the following 
manner. Suppose u satisfies the equation of motion 

M— + C^+K« = Fe‘“" (13) 

dt 2 * 

where M,C and K are matrices derived from finite element analyses of the 
blade and F is a force acting on the blade (it could correspond to a unit 
excitation at the jth node, for example). Let <|>j be the modal displacement 

vector associated with the jth undamped mode, i.e. let <J>j satisfy 

K* . -a? M 4 j O 4 ) 

where coj is the jth resonant frequency of the blade. Then from linear algebra 

we know that the modal displacement vectors are orthogonal with respect to 
M, i.e. 


<j> ^ M 4 j =0 when n = j 

Let C = a M + p K. Then if the response u is 


u = Ue 


i co 1 


(15) 


(16) 


it may be shown (using (14) and (15)) that the coefficients U are given by 


u _y »_l£ 

j kj + icoCj-mjCD 2 


( 17 ) 
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where 


T 

rn j =<|> j M 9 j 

(18a) 

k r m J “i 

(18b) 

and 


c = a m . + pk. 

(18c) 

Thus, the receptances of equation (7), can be found from (17) provided 

we let F be a unit excitation at the jth node and calculate the components of 
U that correspond to the ith node. 


2.3 TUNED SYSTEM THEORY 

The goal of this section is to rewrite the equation that governs blade motion, 
equation (7), in terms of the relative motion between blades. This does two 
things for us. The first is that it becomes clear that for a tuned bladed disk 
all blades respond with the same amplitude (phases differ by a fixed amount) 
and, consequently, the problem can be formulated in terms of a single blade 
or in terms of the motion between a pair of blades. In this section we will 
derive the governing equations in terms of the relative motion between 
blades, since it is these motions that act on the friction dampers and that 
will be related in a later section to the damper forces, fj. 

The tuned system assumption is that all the blades are identical and 
consequently have exactly the same frequencies of vibration (thus, they are 
tuned ). It is also assumed that the blades see exactly the same excitation 
except for the phase of the excitation which differs by a fixed amount, <p, 
from blade to blade, where 


9 = 


"b 


(19) 


and where n nd equals the number of nodal diameters of the response (also 
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equals the engine order of the excitation) and where n^ is the number of 
blades on the disk. 

Consider Fig. 2, a schematic of a section of the bladed disk assembly. In this 
schematic we are looking down on a typical blade in the center and two 
blades on either side. The numbered nodal points are the damper contact 
points. Each nodal point has displacements and forces in the x and y 
directions as indicated in Fig. 3. Because of the phase difference of (19) the 
displacements and forces in the neighboring blades are (quantities with (') 
refer to the blade on the right, those with (") refer to the blade on the left 
and quantities with no (’) refer to the center blade) 

uj = U.e iq> (20a) 

F‘ = Re* (20b) 


Thus, equation (17) can be written in terms of the motions that act on the 
damper between the center and right hand blades, i.e in terms of U 1t U 2 , U 3 ', 

and U 4 ' and in terms of the forces that act on the damper, , F 2 , F 3 \ and 
F 4 \ From (17) and (20) we can write 


Ui +R U f i 


+ R i 2 f 2 




( 21 ) 


We can use (21) for i ■ 1,2. From (20) (a) and (21) 



^[f 0 r r +R y f l +R 12 f 23 +R i3^ +R l 4 f 4 


( 22 ) 


which we can use when i equals 3 or 4. Thus, at this point in the 
formulation we have four equations and eight unknowns (4 displacements and 
4 forces that act on the damper). The goal of the next sections will be to 
relate the damper forces and displacements and to complete the formulation. 


2.4 DAMPER THEORY 


The objective of this section is to simplify the formulation and reduce the 
number of unknowns. The blades and a simplified damper configuration are 
depicted in Fig. 4. Note that in order to simplify the formulation it is 
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assumed that dampers act independently on the upper and lower pairs of 
contact points. It is also assumed that the dampers' inertia are negligible 
(i.e. they behave as springs). Thus, from equilibrium there are two 
independent force vectors since 



(23) 


The damper force can be related to the relative motion across the damper, 
thus the independent displacement quantities of interest are 

w, = Uj - u; (24a) 

W 2 = U 2" U 3 (24b) 

from (21) and (22) we can write 

w i - + tR u + R«-e i ’R 1 4-e i9 R 4 i]<i 

+ [R u+ Ve- l ’R c -' i ' P R4 2 l f 2 (25) 

and 

+ [^22 + ^ 33 ’ C ^ 23 " e < P ^ 32 ^ ( 26 ) 


or 


W = [Wj , w 2 ] T = f 0 r* E + R*F* 


(27) 


where 


F* = [f r f 2 ] T 


and r* E and R* are defined to be consistent with (25) and (26). 


2.5 DAMPER FORCES: P - EM FOR TWO DIMENSIONAL MOTIONS 
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2.5.1 Differences Between One and Two Dimensional Motion: A Special Case 
of Circular Motion 

Consider the case of damper element depicted in Fig. 5 (a) in which the 
motion is only in one direction. The force-displacement relationship is 
depicted in Fig. 5 (b). Note that the damper slips during part of the cycle no 
matter how large the displacement amplitude becomes. Now consider the 
case of a damper element in which the input displacement w is a vector that 
can move in two directions. If the input motion (i.e. the relative motion 
between the blades) is circular then w can be given as 

T 

w = w Q [ cos cot , sin cot] (28) 

and if the damper spring is isotropic, i.e. the damper stiffness matrix is of 
the form 


f = K z and K= k d I where I is the unit matrix (29) 

then the motion of the damper will be circular, but may lag the input by an 
angle, 5. Mathematically, 

llN _ 

z = [ cos(cot - 6) , sin (cot - 5) ] T (30) 

k. 


where 


5 = cos' 1 


f |lN ' 


V w ° k " J 


where pN is the magnitude of the force required for slip. Note that slip only 
occurs if w 0 > pN/k d and that when slip does occur the damper never sticks. 

This is a fundamental difference between friction contact in one and two 
dimensions, i.e. for large motions if the motion is two dimensional the 
damper slips all of the time. 


2.5.2 Establishing the Principal and Minor Axes of a Vector Following an 
Elliptical Orbit 
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For two dimensional motions, the displacement and friction forces at 
contacting points can be represented as two component vectors. In the use 
of receptance formulation, each component of the vector is assumed to be 
sinusoid with an arbitrary phase. In this section it will be shown that such 
vectors follow an elliptical path and that as a consequence the magnitude of 
the vector has a maximum equal to the length of the principal axis of the 
ellipse and has a minimum equal to the length of the minor axis. This section 
summarizes the mathematics necessary to determine the characteristics of 
the ellipse. This information is important in that it provides a basis for 
establishing whether or not the force in a node is large enough to cause slip. 

Considering a periodically varying vector of the form 


f(t)=F 


COSCDt] 

sincttj 


(31) 


where 




y* 


Figure 6 shows that such vector follows an elliptical path whose principal 
direction is at an angle <|> z with respect to the x direction 


Using a singular value decomposition of the matrix F, one may write 


where 


- Jtj 0 

0 X, 





cos<f> z -sin<J> z 
sin<f> 2 cos<{> z 


(32) 


and 
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It is clear that z 1 is a unit vector pointing in the direction of the principal 
axis of the ellipse, is the length of the principal axis, X 2 is the length of 
the minor axis, and <j> s is the phase lag of the motion. 

2.5.3 Elliptical Motion When Slip at Contact Point Does Not Occur 
Consider the case when the contacting point does not slip. Then 

z = w (35) 
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where w is the motion between blades and z is the motion of the damper. 

If w is assumed to be sinusoid, then w can be expressed as the following. 

cosox" 
sincot. 


w xc 

W XX 

Wyc 

w yi 


( 36 ) 


and the friction force becomes 


f(t) = 


k ll k 12 
*22 


W X c w x* 


Wy C Wy, 


coscot 

sincot 


(37) 


Then, using equation (33), it is seen that slip does not occur when X-j<_pN. 


2.5.4 Elliptical Motion at a Fully Slipping Node 

As the periodical motion becomes larger (X-j> pN), slip at the interface 

occurs. The problem of analyzing slip-stick motion for the two-dimensional 
case is very complex. However, if the motion is large enough, one can 
assume that the interface slips all the time. Since w is large, the motion 
across the contact point is given by 


v = w - z r w (38) 

By using equation (33) and rotate the coordinate frame to align with the 
principal axes, v can be expressed as 


v 


a cos0 
b sin0 


where 0 = cot - <f> s and the friction force 

f=-pNv/lv| 


(39) 


becomes 
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a sin8 


0 - 


f(t) = - M-N 


V a 2 sin 2 0 + b 2 cos 2 0 
•f b cos8 

LV a 2 sin 2 0 + b 2 cos 2 0_J 


(40) 


It is seen that f is not a sinusoid. For first order approximation, f may be 
expanded in a Fourier series to obtain its fundamental harmonic components. 
As a result, f can be expressed as 


f (t) = pN 



COS0 

t 

_sin0_ 


(41) 


where 


zr. 


= _L [ 

a sin0 cos0 ( 

T 1 J 
1 0 

V a 2 sin 2 0 + b 2 cos 2 © 

zr, 

-J_ 1 

a sin0 sin0 ( 

T i J 
1 0 

V a 2 sin 2 0 + b 2 cos 2 0 

zr. 

_ 1 

f - b cos0 cos0 

T. J 
1 0 

V a 2 sin 2 0 + b 2 cos 2 0 

zr. 

_ 1 

f - b cos0 sin 0 

T, , 

J / 2 • 2 r\ 


d6 


d0 


d0 


It may be shown that F x c and F y s are zero and F x s and F y c can be found 
explicitly in terms of elliptic integrals. 

2.5.5 Interpolation Method for Estimating Friction Force 

RI/RD 91-230 
A14 



The problem of analyzing slip-stick motion for the two-dimensional case is 
very complex. Since we know the force displacement relationship for the 
damper when the damper is stuck (equation (37)) and when the interface is 
fully slipping (equation (41)), we can set up an interpolation method for 
estimating the harmonic coefficients of the friction force for intermediate 
values of displacement. 

When the motions are small, the friction force in (37) is pure spring force 
which may be expressed as f s (w). When the motions are large, the friction 

force in (41) is pure damping force which may be expressed as ^(w). Then a 

general formulation of interpolation for estimating the friction force can be 
expressed as 


f = s(^ r X 2 ) f s (w) + D(X r X 2 ) f d (w) (42) 

where and X 2 are singular values of the matrix KW and s(X 1 ,X 2 ) and D 
(X 1 ,X 2 ) should meet the following boundary conditions. 


\< jiN 


'sCXj,^) = 1 

pa 1 ,x 2 )=o 


Xj » p.N 


r S(k v X 2 ) = 0 

'p(k v \ 2 ) = 1 


(43a) 


(43b) 


To the present time, the selection of these two interpolation functions is not 
clear. The interpolation functions used in BLDAMP are adapted from the 
solution for one-dimensional motion. They are 

* 

r i x,<^n 

S(X r X 2 ) = jx [0* . 0.5 sin(20*)] Xj > pN (44) 
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* 

where 0 


D < W - i ( i..fiL) 




(45) 


- cos” 1 [1 - 2jiN/X.,]. 
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3. CODE INPUT AND PROGRAM FLOW CHART 


The input requirements for the code are: 

A. Blade Information: 

njj the number of blades. 

Eq engine order of excitation « (# of nodal diameters). 

N m number of modes (all modes correspond to E 0 nodal diameter disk 
modes). 

coj resonant frequency of each mode, i ■ 1, N m . 

<}>j modal displacements for each mode. 

N r number of reference points 

Sj modal stress vector for reference points. 

mj modal mass for each mode. 

Cj modal damping for each mode. 

B. Damper Properties: 

Kdi » N?* and K® , N® ..... damper stiffness matrices for 2 slip loads i = 1,2 
Pj coefficients of friction for upper and lower dampers, 
a, N* N, - a N* and N 2 = (1-a) N* 

C. Excitation Options . 

1. Simulate a tracking plot, i.e. find stress versus frequency of 

excitation for a given input (refer to Fig. 7(a)). 

2. For a given level of input excitation find the peak stress as a 

function of N* (refer to Fig. 7(b)). 

3. For a given damper normal load find the peak stress with the 

damper in place as a function of what the stress would be 
without the damper (refer to Fig. 7(c)). 

The user will be able to select what type of analysis he would like to run. 
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The code will call a subroutine to calculate the undamped response of the 
blade and establish the blade's receptance to the external excitation. There 
will be a default version of the subroutine which only requires the user to 
identify which mode is of interest and the magnitudes of the peak modal 
stresses that would occur if the damper is not present. Calculations will 
then proceed automatically tracing the response of the mode of interest. The 
code will be written so that the function of this subroutine can be replaced 
with a user written version. Thus, the user will be able to specify his own 
excitation model as required. 

The flow chart for the code is indicated in Fig. 8. 
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FIGURE 1. BLRDE SCHEMATIC 
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FIGURE 2. SCHEMATIC OF BLADE ASSEMBLY 
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Blade 


FIGURE 4. DRMPER CONFIGURATION 
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6.0 APPENDIX 


6.2 BLDAMP REQUIRED INPUT DATA AND CODE LISTING 
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BLDAMP input ba ia 


TITLE - 72 Character Title 


I TYPE 


1 = Frequency Response Plot 

\ I 5 tr ® ss vs - No rmal Force Optimum Curve 
o - Performance Curve 


NBLADE , NEODER , NMODES 
EMASS , WMODAL , DAM 


Number of blades, Engine order 
Number of modes 

Generalized mass, Modal frequency 
modal damping ratio 


of excitation 
, Viscous 


} 


(Repeat NMODES times) 


U( 1’D, U(l,3), ... , 0(1, NMODES) 

U(2,l), 0(2,2), 0(2,3), ... , 0(2, NMODES) 


x direction modal 
displacement of damper 
contact point one 
y direction modal 
displacement of damper 
contact point one 


0(8,1), 0(8,2), 0(8,3), ... , 0(8, NMODES) 


NRS , NC0 


Number of reference stresses, 
each reference stress 


Number of components for 


SR( 1 , 1 ) , SR( 1 , 2 ) , ... 
SR(2, 1) , SR( 2 , 2 ) , ... 


, SR(1, NMODES) - Modal stresses 
, SR(2, NMODES) 


SR(NRS*NC0, 1) , SR(NRS*NC0, 2) , 


SR ( NRS*NC0 , NMODES ) 


FNM11 

Kll, K12, K21 , K22 


Lower normal load for contact point 1 
Damper stiffness matrix 


FNM12 

Kll, K12 , K21 , K22 


Opper normal load for contact point 2 
Damper stiffness matrix 


FNM21 

Kll, K12 , K21 , K22 


Lower normal load for contact point 1 
Damper stiffness matrix 


FNM22 

Kll, K12 , K21 , K22 


Upper normal load for contact point 2 
Damper stiffness matrix 


FNMC11,FC0E11 
FNMC12 , FC0E12 
FNMC21 , FC0E21 
FNMC22 , FC0E22 


Lower normal load, friction coefficient 
upper normal load, friction coefficient 
Lower normal load, friction coefficient 
Upper normal load, friction coefficient 


Fraction of normal load at contact point 1 

RI/RD 91-230 
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ALPHA 



IA1 - Type of input excitation data 
If IA1 = 1 Then, 

PSTRESS(l), SVRM(l) - Max undamped stress, modal stress (mode 1) 
PSTRESS(2), SVRM ( 2 ) 


PSTRESS ( NMODES ) , SVRM ( NMODES ) 

If IA1 = 2 Then, 

GF(1) - Generalized force for mode 1 
GF(2) - Generalized force for mode 2 


GF( NMODES) 

If ITYPE = 1 : Frequency Response Analysis 
WMIN, WMAX, DW - Minimum excitation frequency, maximum excitation 

frequency , frequency increment 
FN - Damper normal load (lbs) 

If ITYPE = 2 : Optimum N Curve Analysis 
NP, NO, NONMOD - Stress reference point, component number for 

optimum curve, nonlinear mode number 

If ITYPE = 3 : Performance Curve Analysis 
FN - Damper normal load 

NP, NC, NONMOD - Stress reference point, component number for 

performance curve, Nonlinear mode number 
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PROGRAM BLDAMP 


C 

C BLDAMP CALCULATES THE FORCED RESPONSE OF FR I CT I ONALLY 

C CONSTRAINED TURBINE BLADES AND CAN BE USED TO OPTIMIZE 

C THE DESIGN OF FRICTION DAMPERS 

C TAPE7= I NPUT , TAPE6=0UTPUT 

C 

C AUTHOR - CHIA-HSIANG MENQ 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/NUMB/NBLADE , NEODER f NMODES ,NRS, NCO 
COMMON/MAIN/PHASE 

COMMON/ BLADE /U C 8 , 1 0) , SR ( 1 5 , 1 0) ,EMASS( 10) ,WMODAL( 10) , 

1 DAM (10) 

COMMON / DAMPER / SMTX 1 1(2,2), SMTX 12(2 , 2) , SMTX2 1(2,2), SMTX22( 2 , 2) , 

1 SI (2,2) ,S2(2,2) , FCOE 1 1 , FCOE 1 2 , FC0E2 1 , FCOE22 , F NM 1 1 , FNM 1 2 , 

2 FNM2 1 , FNM22 , FNMC 1 1 , F NMC 1 2 , FNMC 2 1 , FNMC 2 2 , COE 1 ,C0E2 
COMMON/NLOAD/FN 1 , FN2 f ALPHA 

COMMON/ CONVE/ CEPS I 

CHARACTER TITLE*72,FINAME*20, FONAME*20 
C 
C 

WRITE ( * . 5 ) 

5 F0RMAT(//,2X, 'ENTER INPUT FILE NAME',//) 

R EAD ( * , 7 ) FINAME 
7 FORMAT ( A20 ) 

WRITE ( * , 8 ) 

0 FORMAT(//, 2X,‘ 'ENTER OUTPUT FILE NAME'.//) 

READ(*, 7) FONAME 
C 
C 

0PEN(7 , FILE=FINAME , STATUS= 'OLD ' ) 

0PEN(6,FI LE= FONAME , STATUS= ' NEW' ) 

C 

C 

C READ ANALYSIS TYPE 

C 

C I TYPE= 1 FREQUENCY RESPONSE PLOT 

C I TYPE= 2 OPTIMUM CURVE (STRESS VS. NORMAL FORCE) 

C I TYPE=3 PERFORMANCE CURVE (DAMPED VS. UNDAMPED STRESSES) 

C 

READ ( 7 , 10) TITLE 
10 FORMAT ( A72 ) 

READ( 7 , * ) I TYPE 
C 

C READ MODAL INFORMATION AND DAMPER PROPLRlltb 

C 

CALL READMD 
C 

C PRINT INPUT PARAMETERS 

C 

CALL PRINT (TITLE) 

C 

I F ( I TYPE .EQ. 1) THEN 
C 

C CALCULATE FREQUENCY RESPONSE PLOT 

C 

CALL 101 
C 

ELSE IF ( I TYPE . EQ . 2) THEN 
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C 

C 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CALCULATE OPTIMUM CURVE 


ELSE IF ( I TYPE . EQ . 3) THEN 

CALCULATE PERFORMANCE CURVE 

CALL 103 
END IF 
END 

SUBROUTINE READMD 

READMD READS BLADE AND DAMPER INFORMATION 

NBLADE - THE NUMBER OF BLADES 

NEODER - ENGINE ODER OF EXCITATION 

NMODES - THE NUMBER OF MODES; UP TO 10 MODES 

EMASS - MODAL MASS 

WMODAL - MODAL FREQ. 

DAM - MODAL DAMPING RATIO 
U - MODAL DISP. AT CONTACT POINTS 

NR S - THE NUMBER OF REFERENCE STRESS POINTS IS MAX) 

NCO - THE NUMBER OF MODAL STRESS COMPONENTS (3 MAX) 

SR ~ MODAL STRESS AT REFERENCE POINTS 

SMTX11 - DAMPER STIFFNESS MATRIX AT CONTACT POINT ft \ 

AT LOWER BOUND OF THE NORMAL LOAD 
SMTX12 - AT THE UPPER BOUND OF THE NORMAL LOAD 
SMTX21 - DAMPER STIFFNESS MATRIX AT CONTACT POINT ft 7 
AT LOWER BOUND OF THE NORMAL LOAD 
SMTX22 - AT THE UPPER BOUND OF THE NORMAL LOAD 
FNM11 - LOWER BOUND OF THE NORMAL LOAD AT POINT U 1 

FNM12 - UPPER BOUND OF THE NORMAL LOAD AT POINT ft 1 

FNM21 - LOWER BOUND OF THE NORMAL LOAD AT POINT ft 2 

FNM22 - UPPER BOUND OF THE NORMAL LOAD AT POINT ft 2 

F COE 1 1 - FRICTION COEFFICIENT AT CONTACT POINT ft 1 

AT LOWER BOUND OF THE NORMAL LOAD 
FCOE12 - A I THE UPPER BOUND OF THE NORMAL LOAD 
FC0E21 - FRICTION COEFFICIENT AT CONTACT POINT ft 2 
AT LOWER BOUND OF THE NORMAL LOAD 
FCOE22 - AT THE UPPER BOUND OF THE NORMAL LOAD 
FNMC11 - LOWER BOUND OF THE NORMAL LOAD AT POINT *1 

FNMC12 - UPPER BOUND OF THE NORMAL LOAD AT POINT ft 1 

FNMC21 - LOWER BOUND OF THE NORMAL LOAD AT POINT *2 

FNMC22 - UPPER BOUND OF THE NORMAL LOAD AT POINT tt2 

ALPHA = FRACTION OF NORMAL LOAD AT CONTACT POINT 1 
CEPSI = MAX ALLOWABLE ERROR FOR CONVERGENCE 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /NUMB/ NBLADE , NEODER , NMODES, NR S , NCO 
COMMON/MAIN/ PHASE 

COMMON /BLADE/U(8 , 1 0 ) , SR ( 1 5 , 1 0 ) , EMASS ( 1 0) , WMODAl ( JO). 

1 DAM (10) 

COMMON / DAMPE R / SMTX 1 1(2,2) , SMTX 1 2 ( 2 , 2 ) , SMTX2 1 (2,2) , SMTX22(2,2) , 

1 SI (2, 2), S2 (2, 2) , FCOE 1 1 , FCOE 1 2 , FC0E2 1 , FCOE 22 , F NM 1 1 ,FNM12, 

2 FNM2 1 t FNM2 2 t FNMC 1 1 , FNMC 1 2 , FNMC2 1 , FNMC 22 , COE 1 ,C0E2 
COMMON/ NL0AD/FN1 , FN2 .ALPHA 

COMMON /CONVE/ CEP SI 
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c 

c 


c 

c 

c 

c 


CEPS 1 = 1 . 0D-09 

READ (7 , *) NBLADE t NEODER , NMODES 
PHASE = 2.0*3. 1 4 1 59 27 * NEODER/NBLADE 
DO 10 1=1 , NMODES 

10 READ ( 7 , * ) EMASS ( I ) , WMODAL ( I ) , DAM ( I ) 

DO 20 1=1,8 

20 READ ( 7 , * ) (U(I ,J) ,J=1 .NMODES) 

READ ( 7 , * ) NRS.NCO 
NN=NRS*NCO 
DO 40 1=1 ,NN 

40 READ( 7 , * ) ( SR ( I , J) ,J=1 . NMODES ) 

READ( 7 , * ) FNM11 

READ ( 7 f * ) SMTX1 1(1,1 ) , SMTX 1 1 ( 1 , 2 ) , SMTX 11 ( 2 , 1 ) f SMTX11 (2,2) 

READ( 7 , * ) FNM12 

READ (7 , *) SMTX 1 2 ( 1,1), SMTX 1 2 ( 1,2) .SMTX 12(2 , 1 ) .SMTX 12(2, 2) 

READ ( 7 , * ) FNM21 

READ ( 7 , *) SMTX2 1(1,1), SMTX 21 (1.2). SMTX 21 (2,1), SMTX2 1(2.2) 

READ ( 7 , * ) FNM22 

READ ( 7 , * ) SMTX22 ( 1 . 1 ) , SMTX22( 1,2), SMTX22(2, 1 ) , SMTX22(2, 2) 

READ ( 7 , * ) FNMC 1 1 , FCOE 1 1 
READ ( 7, *) FNMC 1 2 , FCOE 1 2 
READ (7,*) FNMC21 .FCOE21 
READ ( 7 , * ) FNMC22 , FC0E22 
READ (7,*) ALPHA 
RETURN 
END 

SUBROUTINE PR I NT ( T I TLE ) 

PRINTS PROBLEM INPUT DATA 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ NUMB/ NBLADE , NEODER , NMODES , NRS , NCO 
COMMON /MAIN/ PHASE 

COMMON/ BLADE /U (8 , 10), SR (15,10), EMASS ( 1 0 ) .WMODAL { 10), 

1 DAM (10) 

COMMON /DAMPER /SMTX 11(2,2), SMTX 12(2,2), SMTX2 1(2,2), SMTX22( 2 , 2) , 

1 SI (2,2) ,S2(2,2) , FCOE 1 1 , FCOE 1 2 , FCOE 2 1 , FCOE22 , FNM 1 1 ,FNM12, 

2 FNM2 1 , FNM2 2 , FNMC 1 1 , FNMC 1 2 , FNMC 2 1 , FNMC22 , COE 1 .C0E2 
COMMON /NLO AD/ FN 1 ,FN2, ALPHA 
COMMON/CONVE/CEPSI 

CHARACTER MD* 4 , TI TLE * 72 
DATA MD/ 'MODE' / 

WRITE(6, 10) 

10 FORMAT (/, 3X, ' — BLDAMP VI. 1 DECEMBER 5,1987 -') 

WRITE(6, 20) 

20 FORMAT(/ ,3X, 'WRITTEN BY C-H. MENQ , GRIFFIN CONSULT I NG , 1 98 7 ' , 

1 / ,3X, 'WARNINGS: 1. THESE CALCULATIONS ARE BASED ON', 

2 / , 1 5X , ' APPROXIMATIONS THAT MAY LEAD TO SIGNIFICANT', 

3 / , 1 5X , ' ERRORS IN CERTAIN CASES — REFER TO USER MANUAL.',/, 

4 / , 1 3X , ' 2 . THE SOLUTION ALGORITHM IS ITERATIVE & MAY NOT', 

5 /, 16X, 'CONVERGE. THIS MAY LEAD TO EXCESSIVE COMPUTER USAGE.') 
WR I TE ( 6 , 30 ) 

30 FORMAT( 1H1 , 2X ,' ******* PROBLEM DESCRIPTION *****♦*') 

WR I TE ( 6 , 40 ) TITLE, NBLADE 

40 FORMAT (/, 2X , A7 2 ,//, 3X , 'NUMBER OF BLADES ' , 26X , ' = ' , I 4 ) 

WR I TE ( 6 , 50 ) NEODER 
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50 FORMAT(3X, 'ENGINE ORDER OF EXC I TATI ON ' . 1 6X . ' = ' , I 4 ) 

WR I TE ( 6 , 60 ) NMODES 

60 FORMAT(3X, 'NUMBER OF MODES ' . 27X . ' = ' . I 4 ) 

WR I TE ( 6 , 90 ) NRS 

90 FORMAT(3X, 'NUMBER OF STRESS REFERENCE POINTS ' , 9X , ' = ' , I 4 ) 

WRITE(6, 1 10) 

110 FORMAT (///,3X, '******* MODAL INPUT DATA ♦***«**',//, 5X ,' MODE ' , 

1 5X, ' FREQUENCY (HZ) 5X ,' GENERALI ZED MASS 5X ,' DAMPING RATIO './ ) 
WRITE (6, 120) (I , WMODAL( I ) , EMASS ( I ) , DAM( I ) , 1=1 .NMODES) 

120 FORMAT ( / , 5X , 1 3 . 6X , 1 PE 1 2 . 5 , 8X , 1 PE 1 2 . 5 . 7X , OPF 1 0 . 4 ) 

WRITE(6. 130) 

130 FORMATC// ,3X, 'MODAL DISPLACEMENTS AT DAMPER CONTACT POINTS',/) 

WR ITE(6, 140) (MD. 1=1 .NMODES) 

140 FORMAT (3X, 'CONTACT' , 6X , B ( A4 . 1 1 X ) ) 

WRITE ( 6 , 145 ) (I , 1=1 .NMODES) 

145 FORMAT (4X. 'POINT' , 7X , 0 ( 1 3 , 1 2 X ) ) 

DO 170 1=1,4 

K=2*I-1 

KK=2*I 

WRITE (6. 150) I , (U( K , J ) , J= 1 .NMODES) 

150 FORMAT (/ , 4X , 1 3 , 4X , 8 ( 1 PE 1 2 . 5 , 3X ) ) 

WRITE (6, 160) (U(KK,J) , J = 1 .NMODES) 

160 FORMAT (11X,8(1PE12.5,3X)) 

170 CONTINUE 

I F ( NR S ,EQ. 0) GO TO 235 
WR I TE ( 6 , 180) 

180 FORMAT(// .3X, 'MODAL STRESSES AT STRESS REFERENCE POINTS'./) 

WR I TE ( 6 , 1 90 ) (MD. 1=1 .NMODES) 

190 FORMAT (2X, REFERENCE' , 5X , 0 ( A4 . 1 1 X ) ) 

WR I TE ( 6 , 1 45 ) (I , 1 = 1 .NMODES) 

NN=NCO*NRS 
DO 230 1 = 1 ,NN 

230 WR I TE ( 6 , 1 50 ) I . ( SR ( I , J ) . J= 1 . NMODES ) 

235 CONTINUE 

WRITE(6. 300) 

300 FORMAT ( 1H1 , IX, '*♦***♦* DAMPER PROPERTIES **♦*♦♦*',//) 

WR ITE ( 6 , 3 1 0 ) 

310 FORMAT(/ ,3X, 'DAMPER CONTACT POINT 1:') 

WR ITE(6,320) FNM11 

320 FORMAT (/, 5X .' LOWER NORMAL LOAD ='.F14.4) 

WR I TE( 6 . 330 ) 

330 FORMAT (/. 5X .' DAMPER STIFFNESS MATRIX:') 

WRITE(6, 340)SMTX1 1(1,1), SMTX 11(1,2) 

WR ITE(6,350) SMTX 1 1(2,1). SMTX1 1(2,2) 

340 FORMAT ( 7X. 'KXX.KXY = ' , 2 ( 2X , F 1 4 . 4 ) ) 

350 FORMAT (7X, 'KYX.KYY = , 2 ( 2X , F 1 4 . 4 ) ) 

WR I TE ( 6 , 360 ) FNM12 

360 FORMAT(/ , 5X ,' UPPER NORMAL LOAD ='.F14.4) 

WRITE(6,330) 

WRITE(6, 340) SMTX 1 2 ( 1,1) ,SMTX12( 1,2) 

WRITE(6, 350) SMTX 12(2, 1 ) .SMTX 12(2, 2) 

WRITE (6, 370) 

370 FORMAT( //, 3X ,' DAMPER CONTACT POINT 2:') 

WR ITE(6,320) FNM21 
WR I TE (6,330) 

WR ITE(6,340) SMTX 2 1 (1,1), SMTX21 (1.2) 

WRITE (6 . 350) SMTX2 1(2,1), SMTX 2 I (2.2) 

WR I TE ( 6 , 360 ) FNM22 
WR I TE ( 6 , 330 ) 

WRITE (6, 340) SMTX22( 1.1) ,SMTX22( 1.2) 
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c 

c 

c 

c 


c 


c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 


WRITE (6, 350) SMTX22(2, 1 ) ,$MTX22(2,2) 

WR I TE (6 , 380 ) 

380 FORMAK// . 10X, COEFFICIENT OF FRICTION 
WR ITE(6,310) 

WR ITE(6,320) FNMC1 1 
WR I TE ( 6 , 390 ) FCOE11 

390 FORMAT ( / # 5X , 'COEFFICIENT OF FRICTION ='.F12.5) 
WRITE(6.360) FNMC12 
WR I TE ( 6 , 390 ) FC0E12 
WR I TE (6,370) 

WR I TE (6,3 20 ) FNMC21 
WR I TE ( 6 , 390 ) FCOE21 
WR ITE ( 6 , 360 ) FNMC22 
WR I TE ( 6 , 390 ) FCOE22 
WRI TE ( 6 , 400) ALPHA 
400 F0RMAT(/,5X, "ALPHA = ' ,1X # F5.2) 

WR I TE ( 6 , 4 1 0 ) CEPSI 
410 FORMAT (/,5X, 'CEPSI = '.1PE12.5) 

RETURN 

END 

SUBROUTINE 101 


101 - PERFORM TRACKING PLOT 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ EXC I T/GF ( 10) .MODE 

WRITE(6 t 10) 

10 FORMAT ( 1H1 , 3X . 'ANALYSIS TYPE SELECTED; FREQUENCY RESPONSE PLOT') 

READ GENERALIZED FORCES FROM STRESS? 

I A 1 = 1 YES 
I A 1 =2 NO 

READ ( 7 , * ) I A 1 
I F ( I A 1 .EQ. 1) THEN 
CALL READS 
ELSE 

CALL READGF 
END IF 


READ FREQUENCY RANGE 


READ( 7 , * ) WMIN.WMAX.DW 
WR I TE ( 6 , 15) 

15 FORMAT ( / , 4X f ' FREQUENCY RANGE;'./) 

WR ITE(6, 20) WMIN.WMAX.DW 
20 F0RMAT(7X, 'MINIMUM FREQUENCY =',F10.3./, 

1 7X. 'MAXIMUM FREQUENCY =* ( F10.3,/. 

2 7X , ' FREQ . INCREMENT =',F10.3,/) 

READ TOTAL DAMPER LOAD 


READ ( 7 , * ) FN 
WR ITE(6,30) FN 

30 F0RMATC4X, 'TOTAL DAMPER LOAD =',F10.3) 
CALL FREQRP (WMI N , WMAX , DW , FN ) 

RETURN 

END 
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CALL READS 
ELSE 

CALL READGF 
ENDIF 
C 

C READ DAMPER NORMAL LOAD 

C 

R EAD ( 7 , * ) FN 
WRITE(6 t 30) FN 

30 FORMAT( /, 7X NORMAL LOAD =‘ .F12.7) 

C 

C READ REFERENCE STRESS, COMPONENT NUMBER, AND MODE 

C 

READ ( 7 , * ) NP , NC , NONMOD 
WR I TE ( 6 , 40 ) NP , NC , NONMOD 
40 FORMAT (/ ,3X, 'PEAK RESPONSE IS BASED ON’,//, 
t 5X , ' STRESS NUMBER =',I3,/, 

2 5X, ' COMPONENT NUMBER =',I3,/, 

2 5X, ' MODE NUMBER =' . 13) 

C 

CALL PER FOR ( NONMOD , FN ) 

RETURN 

END 

C 

SUBROUTINE READGF 
C 

C READS GENERALIZED FORCES 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ EXCIT/GF( 10) .MODE 

COMMON/NUMB/NBLADE .NEODER , NMODES ,NRS , NCO 
DO 20 1=1 .NMODES 
WRITE(6, 10) I 

10 FORMAT ( / , 7X , ’EXCITATION OF MODE NUMBER ' , 1 X , I 2 ) 

READ ( 7 , * ) GF ( I ) 

WR I TE ( 6 , 15) GF ( I ) 

15 FORMAT(/,7X, 'GENERALIZED FORCE = '.1PE15.5,/) 

20 CONTINUE 
RETURN 
END 
C 

SUBROUTINE READS 
C 

C READ PEAK STRESS & MODAL STRESS THEN CALCULATE 

C GENERALIZED FORCE 

C PSTRESS = MAXIMUM UNDAMPED STRESS 

C SVRM = MODAL STRESS 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ EX C I T/GF ( 10) , MODE 

COMMON/NUMB/NBLADE , NEODER . NMODES , NRS , NCO 

COMMON/ BLADE / U ( 8 , 10), SR ( 1 5 p 103 ,EMASS( 10) ,WMODAL( 10), 

1 DAM (10) 

DO 50 1=1 , NMODES 
READ ( 7 , * ) PSTRESS, SVRM 
WR I TE ( 6 , 10) I 

10 F0RMAT(/,7X, 'EXCITATION OF MODE *',1X,I2) 

WRITE (6, 20) PSTRESS, SVRM 

20 FORMATC /, 7X, 'MAXIMUM UNDAMPED STRESS = ' , 1 PE 1 5 . 5 , / . 

1 7X, ' MODAL STRESS = ' , 1 PE 1 5 . 5 . / ) 



GDRM = PSTRESS/SVRM 

GF ( I ) = 2 . 0 * DAM ( I ) *WMODAL ( I ) **2*EMASS ( I ) *GDRM 
50 CONTINUE 
RETURN 
END 
C 

SUBROUTINE FREQRP (WMI N , WMAX , DW , FN ) 

C 

C FIND STRESS VERSUS FREQUENCY OF EXCITATION FOR A GIVEN 

C INPUT. 

C 

C WMIN - SMALLEST FREQUENCY 

C WMAX - LARGEST FREQUENCY 

C DW - REFERENCE FREQUENCY INCREMENT 

C FN - NORMAL LOAD 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / BLADE /U( 8, 1 0 ) , SR ( 1 5 , 1 0 ) , EMASS( 10) ,WMODAL( 10) . 

1 DAM ( 10) 

COMMON/ NUMB /NB LADE t NEODER , NMODES , NRS , NCO 
COMMON/MAIN/PHASE 

COMMON/R EC E PE/RE 1 (2,2),RE2(2,2),RE3(2 I 2),RE4(2,2) 

COMMON /RECEPF/RF(8 , 8 , 2 ) , CF ( 8 , 2 ) 

COMMON /DAMPER / SMTX 1 1(2,2), SMTX12(2,2) , SMTX21 (2,2), SMTX22(2, 2) , 
1 SI (2,2),S2(2,2) , F COE 1 1 , FCOE 1 2 , FC0E2 1 , FC0E22 , F NM 1 1 , FNM 1 2 , 

1 FNM21 , FNM22 , FNMC 1 1 , FNMC 1 2 , FNMC2 1 ,FNMC22,C0E1 , C0E2 
COMMON / NLQAD/ FN 1 , FN2 , ALPHA 
COMMON /EXC IT /GF( 10) , MODE 

CQMMON/GSOLV/X ( 8 ) ,DX(8) ,XMAX(8) ,XMIN(8) ,DELMX(8) , EPSI (8) , 

1 A ( 8 , 8 ) , B ( 8 , 1 ) 

COMMON /ST I CK/ IS , SFN , SW 
COMMON/NITER/ITER 
DIMENSION PS(3) , SM(3 , 2) 

DIMENSION PSS ( 500 ,5,3) ,WW(500) , SMM ( 500 , 5 , 3 , 2 ) 

DIMENSION XI ( 8 ) 

C 

WRITE(6, 10) 

10 FORMAT (/// ,3X, '*♦*** CALCULATE FREQUENCY RESPONSE ♦****') 

C 

CALL GETRE (WMIN ) 

C 

MODE= 1 
CALL DVALUE 

XI ( 1 )=RE1 (1,1 ) +RE4 ( 1,1) 

X I ( 2 ) =RE 1 ( 1 , 2 ) +RE4 (1,2) 

XI (3)=RE1 (2,1 ) +RE4 ( 2 , 1 ) 

XI ( 4 ) =R E 1 (2 , 2)+RE4( 2 ,2) 

XI (5)=RE2( 1 . 1 )+RE3( 1,1) 

XI ( 6) =RE2 ( 1 , 2 ) +RE3 ( 1 ,2) 

XI (7)=RE2(2 , 1 )+RE3(2 , 1 ) 

XI (8)=RE2(2 . 2)+RE3( 2 ,2) 

CALL DAMPR(FN) 

N0EQN=8 
MI TER=20 
IC = 0 
NUMO = 0 
DW 1 = DW 
DW3=DW*3 
I W- 2 

20 W = WMIN 
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30 IC = IC + 1 
40 CALL GETRE(W) 

CALL GETRF(W) 

DO 50 1=1,8 
50 X ( I ) = XI ( I ) 

CALL GSLV ( NOEQN , X , DX , XMAX . XMI N , DELMX , EPS I .20,0,1 ERR . A , B ) 

I F ( I ERR . EQ . 1 ) THEN 
I F ( NUMO .EQ. 0) THEN 
WM I N = 0 . 9 * WM I N 
GO TO 20 
ENDIF 

DW 1 =DW 1 / I W 
W=W~DW 1 
GO TO 40 
ENDIF 

DO 60 1=1,8 
60 XI ( I ) = X ( I ) 

NUMO = NUMQ+ 1 

WW ( NUMO) =W 

DO 80 J=1 . NRS 

CALL OSTRESS(W, J.PS.SM) 

DO 70 JJ=1 f NCO 

PSS ( NUMO , J , J J ) = PS(JJ) 

SMM ( NUMO , J , J J , 1 ) = SM ( J J , 1 ) 

70 SMM ( NUMO . J , J J , 2 ) = SM(JJ,2) 

80 CONTINUE 
C 

I F ( I C .G£. 10) THEN 

IC = IC-10 

ENDIF 

DW1=DW1 *(3.0-8. 0/3 .0*1 TER /MI TER) 

I F ( DW 1 .GE. DW3 ) THEN 

DW 1 =DW3 

ENDIF 

W = W + DW1 
I F ( W .LE. WMAX) THEN 
GO TO 30 
ENDIF 
C 
C 

C PRINT STRESS VS FREQUENCY 

C 

DO 140 J=1 .NRS 
DO 140 K= 1 , NCO 
WR ITE(6, 1 10) J.K 

110 FORMAT( 1H1 , IX, 'FREQUENCY RESPONSE CURVE FOR STRESS NUMBER', 

1 13.' COMPONENT NUMBER I 3 ,//. 6X FREQUENCY 1 OX ,' STRESS ',/ ) 

DO 120 1=1, NUMO 

120 WR I TE ( 6 , 130) WW ( I ) , PSS ( I , J , K ) 

130 FORMAT ( 1 X , 1 PE 1 5 . 5 , 3X , 1 PE 1 5 . 5 ) 

140 CONTINUE 
500 WR I TE ( 6 , 5 1 0 ) 

510 FORMAT (/ ,3X, '**♦** END OF OUTPUT DATA *****•) 

RETURN 

END 

C 

SUBROUTINE OPT I M ( NONMOD ) 

C 

C CALCULATES OPTIMUM NORMAL FORCE CURVE 
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C FN = NORMAL LOAD 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / BLADE / U ( 8 , 10) , SR (15, 10) ,EMASS( 10) , WMODAl. ( 10) , 

1 DAM (10) 

COMMON/ NUMB /NBLADE , NEODER , NMODES , NR S , NCO 
COMMON/MAIN/ PHASE 

COMMON/DAMPER /SMTX1 1(2,2) , SMTX 1 2 ( 2 , 2 ) , SMTX21 (2,2) ,SMTX22(2,2) , 

1 S1(2,2),S2(2,2) , FCOE 1 1 , FCOE 1 2 , FC0E2 1 , FC0E2 2 , FNM 1 1 ,FNM12. 

2 FNM2 1 , FNM2 2 , FNMC 1 1 , FNMC 1 2 , FNMC 2 1 , FNMC 2 2 , COE 1 ,C0E2 
COMMON /NLOAD/FN 1 , FN2 , ALPHA 

COMMON/ RECE PE /R El (2,2), RE2( 2 , 2 ) , RE3 ( 2 . 2 ) , RE4 ( 2 . 2 ) 

COMMON/ R ECEPF /RF (0,8,2) ,CF(8,2) 

COMMON/ FFORCE/ F 1 (2,2),F2(2,2),F3(2,2),F4(2,2) 

COMMON/ EXCIT/GF( 10) .MODE 

C0MM0N/GS0LV1 /XI Cl). DX 1(1) ,XMAX1 ( 1 ) , XMIN 1 ( 1 ) , DELMX 1 ( 1 ) , EPS I 1 ( 1 ) , 
1 A1 ( 1 , 1 ) ,B1 ( 1 . 1 ) 

COMMON /G SO LV/X ( 0 ) ,DX(8) ,XMAX(8) .XMIN (8) . DELMX ( 8 ) .EPS I (8) . 

1 A ( 8 , 0 ) ,B(8, 1) 

COMMON/ ERROR /IR.XI (8) 

COMMON /STICK /I S, SFN, SW 
COMMON/ SW I TCH/ I M 
DIMENSION PS (3) ,SM(3,2) 

DIMENSION PSS( 1 00 , 5 , 3 ) , WW ( 1 00 ) ,SMM( 100,5,3,2) ,FNO( 100) 

DIMENSION SVALUE ( 2 ) , R ANGLE ( 2 ) , PANGLE ( 2 ) 

C 

WR I TE ( 6 , 10) 

10 FORMAT ( / / / , 3X , ****** CALCULATE PEAK RESPONSE *+♦♦*') 

MODE=NONMOD 

XII =WMODAL ( MODE ) ♦ ( 1 . 0-2 . 0 * DAM ( MODE ) ** 2 ) * *0 . 5 
XI I I =0 . 9*XI I 
CALL GETRE1 (XII) 

X( 1 ) =R E 1 (1,1 ) + R E4 ( 1,1) 

X( 2)=RE l ( 1 , 2)+RE4( 1,2) 

X( 3 ) = RE 1 (2.1 ) + R E4 ( 2 , 1 ) 

X ( 4 ) =R E 1 ( 2 , 2) +RE4(2 , 2) 

X ( 5 ) =RE2 ( 1 , 1 ) +RE3 ( 1,1) 

X(6)=RE2( 1 , 2 ) +RE3 ( 1,2) 

X(7)=RE2(2, 1 ) + R E3 ( 2 , 1 ) 

X(8)=RE2(2,2)+RE3(2, 2) 

CALL DVALUE 
CALL DVALUE1 
DFNN=ABS(GF( MODE ) ) 

MI TER=20 
NF = 0 
IC = 0 
NUM0 = 0 
FN = Q . 0 
DFN=DFNN 
20 10=10+1 

NUM0=NUM0+ 1 
I T = 0 
I N = 3 

30 CALL DAMPR(FN) 

XI ( 1 ) = X I I 
DO 40 1=1,8 
40 XI ( I )=X( I ) 

I M= 1 

CALL GSLV1 ( 1 ,X1 , DX 1 ,XMAX1 , XM I N 1 , DELMX 1 ,EPSI 1 .MITER, ITER.O, I ERR. 

1 A 1 , B 1 ) 
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I P ( IERR.EQ. 1 .OR . 1R . EQ. 1 .OR .X 1 ( 1 ) . LE . XI I *0. 99 . OH . 

1 XI ( 1 ) .GT. X I I ♦ 1 . 1 ) THEN 
I T = I T+ 1 

IF(IT .GE. 10) GO TO 100 

FN=FN~DFN*( I N — 1 )/IN 

DFN=DFN/IN 

GO TO 30 

ENDIF 

XII=X1 ( 1 ) 

DO 50 1=1 ,8 
50 X(I)=XI(I) 

I M=2 

CALL GSLV1 ( 1 ,X1 , DX 1 , XMAX1 . XMIN1 , DELMX1 . EPS I 1 . M I TER . I TER , 0 , I ERR , 
1 A 1 , B 1 ) 

W = XI ( 1 ) 

FNO ( NUMO ) = FN 

WW ( NUMO ) =W 

DO 70 J=1 , NRS 

CALL OSTRESS(W, J , PS, SM) 

DO 60 JJ=1 , NCO 
PSSCNUMO, J , JJ)=PS( JJ ) 

SMM ( NUMO , J , J J , 1)=SM(JJ , 1 ) 

60 SMM ( NUMO ,J,JJ,2)=SM(JJ,2) 

70 CONTINUE 

IF(IC.GE. 5) THEN 

IC = I C-5 

ENDIF 

DFN=DFN* (3 . 0-8 . 0/3 . 0*ITER/MITER) 

FN = FN +DFN 
CR I T=W/ X I 1 1 - 1 .0 

I F ( ABS ( CR IT) .GE. 1.0D-07) THEN 
XIII =X 1(1) 

GO TO 20 
ENDIF 
C 

DO 80 1=1,2 
DO 80 J=1 ,2 
FI (I , J) = CF ( I , J ) 

80 F2 ( I , J) = CF(I+2, J) 

CALL ROT ATE ( F 1 , S VALUE , RANGLE , PANGLE ) 

FN 1 = SVALUE ( 1 ) /COE 1 

CALL ROTATE (F2 , S VALUE , R A NGLE , PANGLE) 

FN2 = SVALUE ( 1 )/C0E2 
C 

IF(ALPHA.EQ.Q.O) THEN 
SFN = FN2 
GO TO 90 

ELSE I F( ALPHA . EQ. 1 . 0 ) THEN 
SFN = FN 1 
GO TO 90 
ENDIF 

FN 1 = FN 1 /ALPHA 

FN2 = FN2/( 1 .O-ALPHA) 

I F ( FN 1 . GT . FN2 ) THEN 
SFN = FN 1 
GO TO 90 
ENDIF 
C 

SFN = FN2 

90 FNO ( NUMO- 1 ) =SFN 
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c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


FNO ( NUMO ) = 1 . 1 *SFN 
GO TO 120 
100 WR I TE ( 6 , 110) 

110 FORMAT ( / , 1 X , ' CONVERGEMENT FAILS; STOP COMPUTING | ' ) 

120 CONTINUE 

PRINT PEAK FREQUENCY 

WR I TE ( 6 , 140) 

140 FORMAT ( 1 H 1 , 2X, ' FREQUENCY OF PEAK RESPONSE AS A FUNCTION OF 
1 ' NORMAL FORCE' ,//,4X, 'NORMAL FORCE 5X FREQUENCY ' , / ) 

DO 150 1=1 , NUMO 
150 WRITE(6, 160) FNO(I) f WW(I) 

160 FORMAT (IX, 1P2E15.5) 

PRINT OPTIMUM CURVE 

DO 240 J=1 , NRS 
DO 240 K=1 ,NCO 
WR ITE ( 6,210) J , K 

210 FORMAT( 1H1 , IX, 'OPTIMUM CURVE FOR STRESS NUMBER', 13, 

1 ' COMPONENT NUMBER I 3 ,//, 4X NORMAL FOR CE ' , 9X . ' STR E S S ' , / ) 

DO 220 1=1 .NUMO 

220 WR I TE ( 6 , 230 ) FNO ( I ) . PSS ( I , J , K ) 

230 FORMAT ( 1 X , 1 PE 1 5 . 5 , 3X , 1 PE 1 5 . 5 ) 

240 CONTINUE 
500 WRITE (6,510) 

510 FORMAT ( / , 3X , '**+** END OF OUTPUT DATA ♦**♦*') 

RETURN 

END 

SUBROUTINE PERFOR ( NONMOD , F NN ) 

SUBROUTINE FINDS THE PEAK STRESS WITH THE DAMPER IN PLACE AS 
A FUNCTION OF WHAT THE STRESS WOULD BE WITHOUT THE DAMPER. 


FNN - GIVEN DAMPER NORMAL LOAD 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /BL A DE/U ( 8 , 1 0 ) , SR ( 1 5 , 1 0 ) ,EMASS( 10 ) ,WMODAL( 10 ) . 

1 DAM ( 10 ) 

COMMON/NUMB/NBLADE , NEODER , NMODES , NRS , NCO 
COMMON /MAIN/ PHA S E 

COMMON/DAMPER/ SMTX 11 ( 2 , 2 ) . SMTX 1 2 ( 2 , 2 ) , SMTX 2 1 ( 2 , 2 ) . SMTX 22 ( 2 . 2 ) , 

1 Si (2,2) ,$2(2,2) , F COE 1 1 , F COE 1 2 , FCOE 2 1 , FCOE 2 2 , F NM 1 1 , FNM12, 

2 FNM 2 1 , FNM 22 , FNMC 1 1 , FNMC 1 2 , FNMC 2 1 , FNMC 22 . COE 1 ,C 0 E 2 
COMMON /NLO AD /FN 1 , FN 2 .ALPHA 

COMMON/ RECEPE/ RE 1 ( 2 , 2 ),RE 2 ( 2 , 2 ),RE 3 ( 2 , 2 ),RE 4 ( 2 , 2 ) 

COMMON / RECEPF /RF ( 8 , 0 , 2 ) , CF(B , 2 ) 

COMMON/ FFORCE/F 1 ( 2 , 2 ),F 2 ( 2 , 2 ),F 3 ( 2 , 2 ),F 4 ( 2 , 2 ) 

COMMON/EXCIT/GF( 10 ) .MODE 

COMMON /G SO LV 1 /X 1 ( 1 ) , DX 1 ( 1 ) , XMAX 1(1), XMI N 1 ( 1 ) , Dt: I MX 1 { 1 ) , F PS II ( 1 ) , 


1 A1 ( 1 , 1 ) ,B1 ( 1 , 1 ) 

COMMON /G SOL V/ X ( 0 ) ,DX(8) 


XMAX ( 0 ) , XMI N ( 8 ) ,DELMX(0) .EPS I (B) , 


1 A ( 0 , 0 ) , B ( 8 , 1) 

COMMON / ERROR/ I R, XI (0) 


COMMON / ST I CK/IS,SFN,SW 


COMMON/ SWITCH/ IM , , 

DIMENSION PDN ( 100,5,3), FNO( 100), SMN( 1 00 , 5 . 3 , 2 ) , WW ( 1 00 ) , WWW ( 1 00 ) 
DIMENSION PDD( 100,5,3) ,PUD( 100,5,3) ,GFF( 100) ,SMM( 100,5,3,2) 
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c 

c 

c 

c 


DIMENSION PS ( 3 ) , SM (3,2) 

DIMENSION SVALUE ( 2 ) , R ANGLE ( 2 ) , PANGLE(2) 

WRITE(6, 10) 

10 FORMAT ( / , 3X , '♦**** CALCULATE PERFORMANCE CURVE ♦****') 

MODE = NONMOD 

XI I =WMODAL (MODE ) * ( 1 . 0 - 2 . 0 * DAM ( MODE ) **2) **0 . 5 
XIII =XI 1+0.9 

CALL GETRE1 (XII ) 

X ( 1 ) =RE 1 ( 1 . 1 ) +RE4 ( 1,1) 

X ( 2 ) =RE 1 ( 1 , 2 ) +R E4 ( 1 ,2) 

X ( 3 ) =RE 1 (2, 1 ) + R E4 ( 2 , 1 ) 

X ( 4 ) =RE 1 (2,2) +RE4 (2,2) 

X(5)=RE2( 1 , 1 ) +RE3 ( 1,1) 

X ( 6 ) =RE2 ( 1 , 2 ) +RE3 ( 1 ,2) 

X ( 7 ) =RE2 ( 2 , 1 )+RE3(2, 1 ) 

X(B)=RE2(2 ,2)+RE3(2, 2) 

CALL DVALUE 
CALL DVALUE 1 
CALL DAMPR(FNN) 

M I TER = 20 
NF = 0 
IC = 0 
NUM0=0 
FN = 0 . 0 

FN 1 = ALPHA * FN 

FN2 = (1.0 - ALPHA) * FN 

DFN=ABS(GF( MODE ) ) 

20 IC=IC+1 

NUMO=NUMO+ 1 
I T = 0 
I N = 3 

30 CONTINUE 
X 1 ( 1 ) =XI I 
DO 40 1=1,8 
40 XI(I)=X(I) 

IM= 1 

CALL GSLV1 ( 1 ,X1 ,DX1 ,XMAX1 ,XMIN1 , DELMX 1 , EPS I 1 , M I T ER , I TER , 0 , I ERR, 
1 A 1 , B 1 ) 

IF ( I ERR .EQ. 1 .OR. IR. EQ. 1 .0R.X1 ( 1 ) .LE.XI I *0.99 .OR . 

1 X 1 ( 1 ) .GT. XII*1.1) THEN 
I T= I T+ 1 

I F ( I T .GE. 10) GO TO 100 

FN=FN-DFN*(IN-1)/IN 

DFN=DFN/IN 

FN 1 =A LPHA * FN 

FN 2= ( 1 . O-ALPHA ) * FN 

GO TO 30 

ENDIF 

XII =X 1(1) 

DO 50 1=1 , B 
50 X(I)=XI(I) 

IM = 2 

CALL GSLV 1 ( 1 , X 1 , DX 1 , XMAX 1 , XMIN 1 , DELMX 1 ,EPSI 1 .MITER. I TER , 0 , I ERR , 
1 A 1 , B 1 ) 

W = XI ( 1 ) 
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FNO(NUMO) =FN 

WW ( NUMO ) =W 

DO 60 J = 1 , NR S 

CALL OSTRESS(W,J , PS, SM) 

DO 60 JJ=1 , NCO 

SMN ( NUMO , J , J J , 1 )=SM( JJ , 1 ) 

SMN ( NUMO ,J,JJ,2)=SM(JJ ,2) 

PDN ( NUMO , J , JJ)=PS(JJ) 

60 CONTINUE 

I F ( I C .GE. 5) THEN 

IC = IC-5 

ENDIF 

DFN=DFN*(3.Q-8.0/3.0*ITER/MITER) 

FN=FN+DFN 

FN 1 = ALPHA * FN 

FN2 = ( 1 . O-ALPHA ) * FN 

CR I T =W/X I 1 1 - 1 .0 

I F ( ABS ( CR IT) .GE. 1.0D-7) THEN 
XIII=X1 ( 1 ) 

GO TO 20 
ENDIF 
C 

C PRINT FREQUENCY OF PEAK RESPONSE VS NORMAL FORCE 

C 

WR I TE (6,70) 

70 FORMAT( 1 H 1 , IX, 'FREQUENCY OF PEAK RESPONSE AS A FUNCTION OF'. 
1 ' NORMAL FORCE ' , / / , 3X , ' NORMAL FORCE 5X FREQUENCY './ ) 

WR I TE ( 6,75) (FNO(I) ,WW(I) , 1=1 ,NUMO) 

75 FORMAT ( IX, 1P2E15.7) 

DO 80 1=1,2 
DO 80 J = 1 ,2 
F 1 ( I , J ) = CF ( I , J ) 

80 F 2 ( I , J) = CFU + 2, J) 

CALL ROTATE ( F 1 , SVALUE , R ANGLE , PANGL E ) 

FN 1 = SVALUE ( 1 ) / COE 1 

CALL ROTATE ( F2 , SVALUE , RANGLE , PANGLE ) 

FN2 = SVALUE ( 1 ) /C0E2 
C 

I F ( ALPHA .EQ. 0.0) THEN 
SEN = FN2 
GO TO 90 

ELSE I F (ALPHA .EQ. 1.0) THEN 
SFN = FN 1 
GO TO 90 
ENDIF 

FN 1 = FN 1 /ALPHA 

FN2 = F N 2 / ( 1 .O-ALPHA) 

I F ( FN 1 .GT. FN2 ) THEN 
SFN = FN 1 
GO TO 90 
ENDIF 
C 

SFN = FN2 

90 FNO ( NUMO- 1 )=SFN 
FNO ( NUMO ) = 1 . 1 * SFN 
C 

GO TO 120 
100 WR I TE ( 6 , 110) 

110 FORMAT (/, 1 X CONVERGEMENT FAILS; STOP COMPUTING!' ) 

120 CONTINUE 
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GFF ( 1 )=0.0 
WWW( 1 ) =WW ( NUMO ) 

DO 130 1 = 1, NR S 
DO 130 J= 1 , NCO 
SMM ( 1 , I , J , 1 )=0 . 0 
SMM( 1 f I , J , 2 ) =0 . 0 
PUD ( 1 , I , J)=0.0 
130 PDD ( 1 , X , J)=0.0 
C 

DO 150 I I =2, NUMO 

I 1 =NUMO- I 1 + 2 

GFF ( II) =FNN/FNO( I 1 ) 

WWW( I I )=WW( I 1 ) 

DO 150 1 = 1 , NR S 
DO 140 J=l , NCO 

SMM (II, I , J . 1 ) =SMN (II , I , J, 1 )*GFF( 1 1 ) 

SMM (II , I , J , 2) =SMN ( 1 1 , I , J , 2)*GFF( I I ) 

PUD ( 1 1 , I , J )=PDN( 1 , I , J) *GFF( 1 1 ) 

140 PDD ( 1 1 , I , J ) =PDN ( I 1 , I , J ) *GFF( I I ) 

150 CONTINUE 
C 

C PRINT PEAK FREQUENCY AS A FUNCTION OF GENERALIZED FORCE 

C 

WRITE(6, 160) 

160 FORMAT( 1H1 , IX, ' FREQUENCY OF PEAK RESPONSE AS A FUNCTION OF 

1 'GENERALIZED FORCE ',//, 3X GENERAL l ZED FORCE , 6X, 

2 'FREQUENCY',/) 

DO 170 1=1, NUMO 

170 WR ITE ( 6 , 180) FNO ( I ) , WWW ( I ) 

180 FORMAT (3X. 1 PE 1 5 . 7 , 4X , 1 PE 1 5 . 7 ) 

C 

C PRINT PERFORMANCE CURVE 

C 

DO 250, J=1 , NR S 
DO 250 K = 1 , NCO 
WR I TE ( 6 , 220 ) J,K 

220 FORMAT( 1H1 , IX, 'DAMPER PERFORMANCE CURVE FOR STRESS NUMBER', 13, 

1 ' COMPONENT NUMBER ' ,I3,//,4X, 'GEN. FORCE/NORMAL FORCE', 1 OX , 

2 ' STRESS * ,/) 

DO 230 1=1 , NUMO 

230 WR I TE ( 6 , 240 ) PUD ( I , J , K ) , PDD ( I , J , K ) 

240 FORMAT (6X, 1 PE 1 5 . 5 , 5X , 1 PE 1 5 . 5 ) 

250 CONTINUE 
500 WR ITE(6,510) 

510 FORMAT (/,3X, '****♦ END OF OUTPUT DATA ***♦♦') 

RETURN 

END 

C 

SUBROUTINE OSTRESS(W,NREFER,PS,SM) 

C 

C CALCULATE STRESS INFORMATION FOR A REFERENCE POINT 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ BLADE/U(8 , 1 0 ) , SR ( 1 5 , 1 0 ) ,EMASS( 10) ,WMODAL( 10) , 

1 DAM (10) 

COMMON / NUMB /NBLADE , NEODER , NMODES , NRS , NCO 

COMMON/MAIN/PHASE 

COMMON /RECEPF/RF (8,8,2) ,CF(8 , 2) 

COMMON/ EX C I T/GF ( 10) , MODE 
DIMENSION SS( 10) ,FS( 10) 
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DIMENSION SM (3,2), PS ( 3 ) 

C 

DO 100 I =1 , NCO 
DO 10 Jl = l , NMODES 
N 1 = NREFER-1 

10 SS( J 1 )=SR(NC0*N1 + 1 , J 1 ) 

C 

CALL RESPNSE(SS,W, SM( I , 1 ) , SM ( I , 2 ) ) 

C 

DO 30 K= 1 , 8 
DO 20 J 2= 1 , NMODES 
20 F S ( J 2 ) = U ( K , J 2 ) 

C 

CALL STRESS (SS.FS. 1 . OD+O , W , SCA , S SA ) 

SM(I.I) = SM (1,1) - CF ( K , 1 ) * SCA + CF(K,2)*SSA 
SM ( I t 2 ) = SM ( I , 2 ) - CF ( K , 2 ) * SCA - CF(K,1)*SSA 
30 CONTINUE 

PS(I)= ( SM ( I , 1)«‘*2-*-SM(I,2)**2)**0.5 
100 CONTINUE 
RETURN 
END 
C 

SUBROUTINE DVALUE 
C 

C SPECIFY NUMERICAL VALUES FOR THE PARAMETERS USED 

C IN GSOLV 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / BLADE / U ( 8 , 1 0 ) , SR ( 1 5 , 1 0 ) ,EMASS( 10) ,WMODAL( 10) , 

1 DAM (10) 

COMMON /GSOLV / X ( 0 ) , DX ( 8 ) , XMAX(0 ) , XMIN (8 ) , DELMX ( 8 ) , EPS I ( 8 ) , 
1 A ( 8 , 8 ) , B ( 8 , 1) 

COMMON /EXCIT/GF( 10) .MODE 
COMMON/MAIN/PHASE 

COMMON/ RECEPE/ RE 1 ( 2 . 2 ) , RE2 ( 2 t 2 ) . RE3 ( 2 . 2 ) . RE4 ( 2 . 2 ) 

COMMON/CONVE/CEPSI 

DIMENSION XX ( 8 ) 

CALL GETRE (WMODAL (MODE ) ) 

XX ( 1 )=RE1 (1,1 ) +RE4 ( 1,1) 

XX ( 2 ) =RE 1 ( 1 , 2 ) +RE4 ( 1 ,2) 

XX ( 3 ) =RE 1 (2 , 1 )+RE4(2. 1 ) 

XX ( 4 ) =RE 1 (2,2) +RE4 (2,2) 

XX ( 5 ) =RE 2 ( 1 , 1 ) +R E 3 ( 1,1) 

XX(6)=RE2( 1 , 2)+RE3( 1 , 2) 

XX ( 7 ) =RE2 ( 2 , 1 )+RE3(2, 1 ) 

XX ( 8 ) =RE2 ( 2 , 2)+RE3( 2 ,2) 

DO 10 1=1,2 
J 1 =4*1 -3 
J2=4*I-2 
J3 = 4 * I - 1 
J4 = 4 * I 

XX2=(XX(J1)**2+XX(J2) + *2+XX(J3)**2+XX(J4)^*2)*+0.5 
EPSI (J1 )=CEPSI*XX2 
EPSI (J2)=EPSI (J 1 ) 

EPSI ( J3)=EPSI ( J 1 ) 

EPSI ( J4) =EPSI ( J 1 ) 

DX( J 1 )=0 . 00000 1 *XX2 
DX ( J 2 ) = DX ( J 1 ) 

DX ( J3 ) =DX ( J 1 ) 

DX ( J4 ) =DX ( J 1 ) 
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XMAX ( J 1 ) = 1 . OD 1 5 
XMAX(J2) =XMAX ( J 1 ) 

XMAX( J3)=XMAX( J 1 ) 

XMAX ( J 4 ) =XMAX ( J 1 ) 

XMI N(J1)=“10000. *XX2 
XMIN( J2)=XMIN( J 1 ) 

XMI N ( J 3 ) =XMI N ( J 1 ) 

XMI N ( J4 ) = XMI N ( J 1 ) 

DELMX ( J 1 ) =0 . 1 *XX2 
DELMX ( J 2 ) =DELMX ( J 1 ) 

DELMX ( J 3 ) =DELMX ( J 1 ) 

10 DELMX ( J4 ) -DELMX ( J 1 ) 

RETURN 

END 

C 

SUBROUTINE DVALUE1 
C 

C SPECIFY NUMERICAL VALUES FOR THE PARAMETERS USED 

C IN GSOLV1 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / BLADE / U ( 8 , 1 0 ) , SR ( 1 5 , 1 0 ) , EMASS( 10) ,WMODAL( 10) , 

1 DAM ( 10) 

COMMON /GSOLV 1 /X 1 ( 1 ) , DX 1 ( 1 ) ,XMAX1 ( 1 ) .XMINI ( 1 ) , DELMX 1 ( 1 ) . EPS 11(1), 
1 A1(1,1),B1(1,1) 

COMMON /EXCIT/GF( 10) .MODE 
EPS I 1 ( 1 ) =0.01 /WMODAL( MODE) 

DX 1 ( 1 )=0.01*DAM( MODE ) +WMODAL ( MODE ) 

XMAX 1 ( 1 )=1 .0015 

XMINI ( 1 )=-10000.0 

DELMX 1 ( 1 )=0. 1 *WMODAL ( MODE ) 

RETURN 

END 

C 

SUBROUTINE DAMPR(FN) 

C 

C CALCULATE DAMPER STIFFNESS AND FRICTION COEFF. 

C AT A GIVEN NORMAL LOAD 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / DAMPER /SMTX1 1(2,2). SMTX1 2( 2 , 2) , SMTX2 1(2,2), SMTX22(2 , 2 ) , 

1 SI (2,2) ,52(2,2) , F COE 1 1 , FCOE 1 2 , FC0E2 1 ,FCOE22,FNM1 1 , FNM 1 2 , 

1 FNM21 , FNM22 , FNMC 1 1 , FNMC 1 2 , FNMC2 1 , FNMC22 , COE 1 ,C0E2 
COMMON /N LOAD/ FN 1 , FN2 . ALPHA 
C 

FN 1 = ALPHA*FN 
FN2 = ( 1 . O-ALPHA ) * FN 
C 

BATA1 = ( FN 1 - FNMC 1 1 ) / ( FNMC 1 2 -FNMC 1 1 ) 

I F ( BATA 1 .GT. 1 . 0 ) THEN 
BATA 1 = 1.0 

ELSE I F ( BATA 1 . LT . 0.0) THEN 
BATA! = 0.0 
END I F 
C 

BATA 2 = ( FN 2- FNMC 2 1 ) / ( FNMC 2 2- FNMC2 1 ) 

I F ( BATA 2 . GT . 1 . 0 ) THEN 
BATA2 = 1.0 

ELSEI F ( BATA2 . LT. 0. 0)THEN 
BATA 2 = 0.0 
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ENDI F 
C 

COE 1 = ( 1 . 0-BATA 1 ) *FCOE 1 1 + BATA1*FCOE12 

COE2 = ( 1 .0-BATA2)*FCOE21 + BATA2*FCOE22 
BATA 1 = (FN1 -FNM11 ) / ( FNM1 2-FNM1 1 ) 

I F ( BATA 1 .GT. 1 . 0 ) THEN 
BATA 1 = 1.0 

ELSEIF (BATA 1 . LT . 0 . 0 ) THEN 
BATA 1 = 0.0 
ENDI F 
C 

BATA2 = ( FN2-FNM2 1 ) / ( FNM2 2- FNM2 1 ) 

I F ( BATA2 . GT . 1 . 0 ) THEN 
BATA 2 = 1.0 

ELSEIF (BATA 2. LT. 0.0) THEN 
BAT A 2 = 0.0 
ENDI F 
C 

SH1.1) = ( 1 . 0-BATA 1 ) * SMTX 1 1(1,1 ) +BATA 1 * SMTXI 2( 1,1) 

S 1 ( 1 , 2 ) = ( 1 . 0-BATA 1 ) *SMTX11 ( 1 , 2)+BATA 1 *SMTX 1 2( 1 , 2) 

SI (2 , 1 ) = ( 1 . 0-BATA 1 )* SMTXI 1(2,1 ) + BATA 1 ♦SMTXI 2 (2 . 1 ) 

SI (2,2) = ( 1 ,0-BATAI ) * SMTX 1 1 (2, 2)+BATA1* SMTX 12(2,2) 

C 

S2( 1 . 1 ) = ( 1 . 0-BATA 2) *SMTX2 1(1,1 ) + BATA 2 * SMTX22 ( 1,1) 

S2( 1 ,2) = (1 . 0-BATA2) * SMTX 2 1 ( 1 , 2 ) + BAT A 2 * SMTX2 2 ( 1,2) 

S2 ( 2 , 1 ) = ( 1 .0-BATA2)*SMTX21 (2, 1 ) +BATA2* SMTX22 ( 2 . 1 ) 

S2 ( 2 , 2 ) = ( 1 ,0-BATA2)*SMTX21 ( 2 , 2 ) +BATA2* SMTX22 ( 2 . 2 ) 

C 

RETURN 

END 

C 

FUNCTION F ( I EQ , NOEQN , X ) 

C 

C A SET OF NONLINEAR EQUATIONS WHICH DESCRIBE THE MOTION 

C AT THE CONTACTING POINT 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /DAMPER /SMTXI 1(2,2), SMTX 12(2,2), SMTX 2 1(2,2), SMTX22( 2 , 2) , 
1 SI (2,2) ,S2(2,2) , FCOE 1 1 , FCOE 1 2 , FC0E2 1 , FCOE22 , FNM 1 1 , FNM 1 2 , 

1 FNM21 , FNM22 , FNMC 1 1 , FNMC 1 2 , FNMC2 1 ,FNMC22,C0E1 ,C0E2 
COMMON/ NLOAD/FN 1 , FN2 .ALPHA 
COMMON / EXCIT/GF( 10) .MODE 

COMMON/ RECEPE/ RE 1 (2,2).RE2(2.2),RE3(2,2),RE4(2,2) 

COMMON/ REC EPF /RF(8,8,2) ,CF (8 , 2) 

COMMON/FFORCE/F 1 ( 2 , 2 ) . F2 ( 2 , 2 ) , F3 ( 2 , 2 ) , F4 ( 2 , 2 ) 

COMMON /MA IN/ PHASE 
COMMON /STICK/IS, SFN.SW 
DIMENSION X(B) ,W1 (2,2) , W2 ( 2 , 2 ) 

DIMENSION CFF(8, 2) ,FF2(2, 2) , FF 1 ( 2 , 2) 

W1 ( 1 , 1 ) = X( 1 ) 

W1 ( 1 , 2) = X ( 2 ) 

W1 (2 , 1 ) = X ( 3 ) 

W1 (2,2) = X ( 4 ) 

W2( 1 , 1 ) = X ( 5 ) 

W2( 1 , 2) = X ( 6 ) 

W2(2, 1 ) = X ( 7 ) 

W2 (2,2) = X ( 8 ) 

C 

IF(IS.EQ. 1 ) THEN 
CALL AB ( S 1 , W 1 , 2 , F 1 ) 
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CALL AB ( S 2 , W2 , 2 , F 2 ) 

GO TO 10 
END I F 

CALL FR I FORC ( W 1 ,S1 ,C0E1 ,FN1 , F 1 ) 

CALL FRIF0RC(W2, S2 , C0E2 ,FN2,F2) 

10 CALL PHASLAG(F1 .PHASE, F4) 

CALL PHASLAG(F2, PHASE, F3) 

C 

DO 20 1=1,2 
DO 20 J= 1 . 2 
CF ( I , J}= F 1 ( I , J ) 

CF ( I + 2 , J ) = F 2 ( I , J) 

CF(I+4, J) = F3 ( I , J) 

20 CF ( 1+6 , J ) = F4(I , J) 

C 

CALL PHASLAG(F1 .-PHASE, FF1 ) 

CALL PHASLAG ( F2 , -PHASE , FF2 ) 

C 

DO 30 1=1,2 
DO 30 J=1 .2 
CFF(I,J)= -FFI(I.J) 

CFF ( I + 2 , J ) = -FF2 ( I , J) 

CFF(I+4, J) = -F 2 ( I , J) 

30 CFF ( 1 + 6 , J ) = -FI (I , J) 

C 

DO 50 1=1 .2 

W1(I, 1 ) =RE 1 ( I . 1)+RE4(I,1) 

W 1 (I , 2 ) =R E 1 (I , 2 } +R E4 ( I ,2) 

DO 40 K=1 ,8 

W1(I,1)= W1(I, 1)-RF(I,K, 1 ) *CF ( K , 1 ) +RF ( I ,K,2)*CF(K,2) 

W1 ( I , 1 ) = W1 ( I , 1 )+RF( 1+6, K , 1 ) *CFF ( K . 1 ) -R F ( I +6 , K , 2 ) *CF F ( K , 2 ) 
W 1 ( I , 2 ) = W1 (I , 2 ) -R F ( I ,K,2)*CF(K, 1 ) -R F ( I ,K, 1 )*CF(K,2) 

40 W 1 ( I , 2 ) = W1(I ,2)+RF(I+6,K,2)*CFF(K, 1 )+RF(I+6,K, 1 ) * CF F ( K , 2 ) 
50 CONTINUE 
C 

DO 70 1=1,2 

W2( I , 1 )= RE2 ( I , 1 ) +RE3 ( 1,1) 

W2 ( I , 2 ) = R E2 ( I ,2)+RE3(I ,2) 

DO 60 K=1 ,0 

W2( I , 1 )= W2( I , 1 )-RF( 1 + 2. K, 1 )*CF(K , 1 )+RF( I +2. K , 2) *CF(K, 2) 
W2( I , 1 )= W2( I , 1 )+RF( 1+4 , K , 1 ) *CFF(K , 1 ) -RF( 1+4, K, 2) *CFF (K, 2) 
W2( I , 2)= W2( I , 2)-RF ( I +2 , K , 2 ) *CF ( K , 1 ) -RF ( I +2 , K , 1 ) *CF ( K . 2 ) 

60 W2 (I , 2 ) = W2(I.2)+RF(I+4.K,2)*CFF(K, 1)+RF(I+4,K, 1)*CFF(K,2) 
70 CONTINUE 
C 

IF(IEQ.EQ.I) THEN 
F=X ( 1 )-W1 (1,1) 

ELSE I F ( I EQ . EQ . 2 ) THEN 
F=X ( 2 ) -W1 (1,2) 

ELSEIF( IEQ. EQ. 3) THEN 
F=X ( 3 ) -W 1 (2,1) 

ELSE I F ( I EQ . EQ . 4 ) THEN 
F=X ( 4 ) -W 1 (2,2) 

ELSEX F ( I EQ . EQ . 5 ) THEN 
F=X(5)-W2( 1,1) 

ELSE I F ( IEQ.EQ.6) THEN 
F=X(6)-W2( 1 ,2) 

ELSE I F ( I EQ . EQ . 7 ) THEN 
F = X ( 7 ) -W2 ( 2 , 1 ) 

ELSE I F ( I EQ . EQ . 8 ) THEN 
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F = X ( 8 ) -W2 ( 2,2) 

ENDIF 

RETURN 

END 

C 

FUNCTION FI ( I EQ , NOEQN , X 1 ) 

C 

C A NONLINEAR EQUATION WHICH DEFINES THE PEAK RESPONSE 

C OF AN OPTIMAL CURVE 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / EXCIT/GF( 10) , MODE 

COMMON/ BLADE /U ( 8 , 1 0 ) , SR ( 1 5 , 1 0 ) , EMASS ( 10) ,WMODAL( 10) , 

1 DAM (10) 

COMMON /RECEPE /RE 1 (2,2) ,RE2(2 # 2) ,RE3(2,2) ,RE4(2.2) 

COMMON/ RECEPF /RF (8,8,2) , CF(0 , 2) 

COMMON /GSOLV/X (8) ,DX(8) ,XMAX(8) ,XMIN(8) ,DELMX(8) .EPS I (8) , 

1 A(B,8) ,B(8,1) 

COMMON/ ERROR / I R , X I (8) 

COMMON/ ST I CK/ I S , SFN , SW 
COMMON/ SWI T CH/ I M 
DIMENSION XI ( 1 ) 

C 

W = XI ( 1 ) 

IT = 30 

I F ( I S . EQ . 1 ) THEN 

IT = 50 

ENDIF 

I F ( I M .EQ. 1) THEN 
CALL GETRE(W) 

ELSE I F ( I M .EQ. 2) THEN 
CALL GETRE 1 ( W) 

ENDIF 

CALL GETRF(W) 

CALL GSLV ( 8 , X I , DX , XMAX , XM I N , DE LMX , E P S I , I T , 0 , i R , A , B ) 

I F ( IR . EQ. I ) GO TO 10 
I F ( I EQ . EQ . 1 ) THEN 

F 1 =F9 ( W , .001 * DAM ( MODE )*WMODAL( MODE) .XI ) 

ENDIF 
10 RETURN 
END 
C 

FUNCTION F9 ( W , DW , X I ) 

C 

C CALCULATE THE SLOPE OF AN OPTIMAL CURVE AT A GIVEN FREQUENCY 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON / BL ADE /U(8 , 10) , SR ( 15 . 10) ,EMASS( 10) .WMODAl ( 10) . 

1 DAM (10) 

COMMON/ NUMB /NBLADE , NEODER , NMODES , NRS , NCO 
COMMON/MAIN/PHASE 
COMMON/ EX C I T / GF ( 10) .MODE 
COMMON/RECEPF/RF (8 , 8,2),CF(8,2) 

COMMON /GSOLV/X ( 8 ) , DX ( 8 ) . XMAX (0) , XMl N( 8 ) , DELMX ( 8 ) .EPSI 18) , 

1 A(8 ,8) ,B(8, 1) 

COMMON /OPT I MUM/ NP , NC 
DIMENSION XI (8) ,WW(2) t FF(2) 

DIMENSION PS (3) , SM( 3 , 2 ) 

C 

WW ( 2 ) = W+0 . 0 1 *DW 
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WW( 1 ) =W~0 . 0 1 *DW 
DO 10 11=1,2 
CALL GETRE ( WW( II)) 

CALL GETRF ( WW ( I I ) ) 

X( 1) = XI ( 1) 

X( 2) = XI (2) 

X ( 3 ) = XI (3) 

X ( 4 ) = XI (4) 

X ( 5 ) = X I ( 5 ) 

X ( 6 ) = XI ( 6 ) 

X ( 7 ) = X I ( 7 ) 

X ( 8 ) = XI ( 8 ) 

CALL GSL V ( 8 , X , DX , XMAX , XM I N , DELMX . EPS I ,50,0,1 ERR , A , B ) 

CALL OSTRESS(WW( II) ,NP,PS,SM) 

10 FF ( 1 1 ) =PS ( NC ) 

F9=(FF(2)-FF( 1 ) )/(0.02*DW) 

RETURN 

END 

C 

SUBROUTINE GETRE (W) 

C 

C CALCULATE VIBRATORY DISPLACEMENTS OF THE FOUR CONTACT 

C POINTS WHEN ASSUMING ZERO NORMAL LOAD (NO FRICTION FORCE) 

C AND MULTI-MODE EXCITATION 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ NUMB / NBLADE .NEODER ,NMODES , NR S.NCO 

COMMON/ BLADE / U ( 8 , 1 0 ) , SR ( 1 5 , 1 0 ) ,EMASS( 10) ,WMODAL( 10) , 

1 DAM (10) 

COMMON /R EC EPE /RE 1(2,2) ,RE2(2,2),RE3(2,2) ,RE4(2,2) 

COMMON/ FFORCE/ FI (2,2),F2(2,2),F3(2,2),F4(2,2) 

COMMON/MA IN/ PHASE 
COMMON /EXCIT/GF( 10) .MODE 
DIMENSION RE5(2,2) ,RE6(2,2) 

DIMENSION GD( 10) 

C 

DO 50 1=1,2 
DO 10 J=1 .NMODES 
10 GD ( J ) = U(I , J) 

CALL RESPNSE(GD,W, RE1 (I,1),RE1(I,2)) 

DO 20 J= 1 , NMODES 
20 GD ( J ) = U( 1+2 , J ) 

CALL RESPNSE ( GD , W , RE2 ( I , 1 ) , RE2 ( I , 2 ) ) 

DO 30 J=1 .NMODES 
30 GD ( J ) = U ( I +4 , J ) 

CALL RESPNSE(GD,W, RE5( I , 1 ) ,RE5(I ,2) ) 

DO 40 J=1 .NMODES 
40 GD ( J ) = U ( I +6 , J ) 

50 CALL RESPNSE ( GD ,W,RE6(I , 1 ) ,RE6(I ,2) ) 

CALL PHASLAG(RE5, -PHASE, RE3) 

CALL PHASLAG(RE6, -PHASE, RE4) 

RETURN 

END 

C 

SUBROUTINE GETRE 1 (W) 

C 

C CALCULATE VIBRATORY DISPLACEMENTS OF THE FOUR CONTACT 

C POINTS WHEN ASSUMING ZERO NORMAL LOAD (NO FRICTION FORCE) 

C AND SINGLE-MODE EXCITATION 

C 
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IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/NUMB/NBLADE , NEODER , NMODES # NRS , NCO 
COMMON/BLADE/U ( 8 , 10) , SR( 15 , 10) , EMASS( 10) .WMODALt 10) , 

1 DAM (10) 

COMMON /RECE PE /RE 1(2,2),RE2(2,2),RE3(2,2).RE4(2,2) 

COMMON/ FFORCE/F 1 (2,2),F2(2,2),F3(2,2),F4(2,2) 
COMMON/MAIN/ PHASE 
COMMON/ EXC IT /GF ( 10) .MODE 
DIMENSION RE5(2,2) , RE6( 2 , 2 ) 

DIMENSION GD( 10) 

C 

DO 50 1=1 .2 
DO 10 J=1 , NMODES 
10 GD ( J ) = 0.0 

GD ( MODE) =U ( I .MODE) 

CALL RE$PNSE(GD,W,RE1 (I . 1 ) , RE 1 ( I . 2 ) ) 

DO 20 J=1 .NMODES 
20 GD ( J ) = 0.0 

GD ( MODE ) =U ( I + 2 . MODE ) 

CALL RESPNSEC GD , W , R E 2 ( I . t ) , RE 2 ( I . 2 ) ) 

DO 30 J=1 .NMODES 
30 GD ( J ) = 0.0 

GD ( MODE ) =U( 1+4 .MODE) 

CALL RESPNSE(GD,W,RE5( I . 1 ) ,RE5( I .2) ) 

DO 40 J=1 .NMODES 
40 GD ( J ) = 0.0 

GD(MODE)=U( 1+6 . MODE ) 

50 CALL RESPNSEC GD . W . RE6 ( I , 1 ) ,RE6(I ,2) ) 

CALL PHASLAGCRE5, -PHASE, RE3) 

CALL PH A SLAG ( RE6 , -PHASE , R E4 ) 

RETURN 

END 

C 

SUBROUTINE GETRF(W) 

C 

C CALCULATE RECEPTANCE MATRIX 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/BLADE/U (8. 10) ,S( 15. 10) ,EMASS( 10). WMODAL (10), 

1 DAM (10) 

COMMON/NUMB/NBLADE , NEODER , NMODES , NRS . NCO 
COMMON/MA I N/ PHASE 

COMMON/ RECEPF/RF (8 , 8 , 2) . CF (8 , 2 ) 

DIMENSION DSHAPEC 10) , FSHAPE( 10) 

C 

DO 20 1 = 1 ,8 
DO 20 J = 1 .8 
C 

DO 10 K= 1 , NMODES 
DSHAPE(K) = U( I .K) 

10 FSHAPE(K) = U ( J • K ) 

C 

CALL RECEPTNC D SHAPE , F SHAPE . 1 .0D + 0,W,RF(I .J, l),fU(I ,J,2)) 
20 CONTINUE 
RETURN 
END 

Q 

SUBROUTINE RECEPTNC D SHAPE , F SHAPE .FC.W.RC.RS) 

C 

C CALCULATE RECEPTANCE (DISPLACEMENT) 
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c 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/ BLADE / U ( 8 ,10),S(15,10), EMASS( 10), WMODAL (10), 

1 DAM (10) 

COMMON/ NUMB /NBLADE , NEODER , NMODES , NR S , NCO 
DIMENSION DSHAPEC 10) , FSHAPE( 10) 

RC = 0.0 
RS = 0.0 

DO 10 J=1 , NMODES 

GFC = F SHAPE ( J ) * FC/ EMASS ( J ) /WMODAL ( J ) * * 2 
DENI = 1.0 - (W/ WMODAL ( J ) ) ♦ ♦ 2 
DEN2 = 2 . 0*DAM( J ) ♦W/ WMODAL ( J ) 

RC = RC + D SHAPE (J) * DEN 1 *GFC / (DENI **2+DEN2**2 ) 

RS = RS + D SHAPE (J ) * (DEN2*GFC) / (DENI **2+DEN2**2) 

10 CONTINUE 
RETURN 
END 
C 

SUBROUTINE STRESS ( SSHAPE , FSHAPE . FC , W. SC . SS ) 

C 

C CALCULATE RECEPTANCE (STRESS) 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON /BLADE / U ( 8 . 10),S(15,10) , EMASS ( 10) .WMODAL ( 10) . * 

1 DAM (10) 

COMMON/NUMB/NBLADE . NEODER , NMODES , NRS , NCO 
DIMENSION SSHAPE( 10) ,FSHAPE( 10) 

C 

SC = 0.0 
SS = 0.0 

DO 10 J=1 .NMODES 

GFC = FSHAPE (J ) * FC/ EMASS (J ) /WMODAL (J ) 

DENI = 1.0 - (W/ WMODAL ( J ) ) * * 2 
DEN2 = 2 . 0*DAM ( J ) *W/WMODA L ( J ) 

SC = SC + SSHAPE (J ) * DEN 1 *GFC / (DENI **2+DEN2**2) 

SS = SS + SSHAPE(J)* (DEN2*GFC)/(DEN1 **2+DEN2**2) 

10 CONTINUE 
RETURN 
END 
C 

SUBROUTINE RE SPNSE ( GD , W , RC , R S ) 

C 

C CALCULATE VIBRATORY RESPONSE OF A NODE POINT WHOSE 

C MODAL DISPLACEMENT/STRESS IS GIVEN 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/NUMB/NBLADE . NEODER . NMODES . NRS , NCO 

COMMON /BLADE /U(8 . 10). SOB, 10) . EMASS ( 10) ,WMODAL( 10) . 

1 DAM (10) 

COMMON/ EXC I T/GF ( 10) .MODE 
DIMENSION GD( 10) 

C 

RC = 0.0 
RS = 0.0 

DO 10 J=1 .NMODES 

GFC = GF ( J ) / EMASS ( J ) /WMODAL ( J ) * * 2 
DENI = 1.0 - (W/WMODAL( J ) ) **2 
DEN 2 = 2 . 0* DAM ( J ) *W/WMODAL ( J ) 

RC = RC + GD( J) *DEN1 *GFC/(DEN1 **2+DEN2**2) 

RS - RS + GD(J)*(DEN2*GFC)/(DEN1**2+DEN2**2) 
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10 CONTINUE 
RETURN 
END 
C 

SUBROUTINE PHASLAG ( F 1 , PHASE , F2 ) 

C 

C CALCULATE PHASE ANGLE BETWEEN ADJACENT BLADES 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-2) 

DIMENSION FI (2,2) ,F2(2,2) 

C 

C = COS ( PHASE ) 

S = S I N ( PHASE ) 

F2( 1 , 1 ) = -C*F 1 ( 1.1) + S*F 1 ( 1 . 2 ) 

F 2 ( 2 . 1 ) = -C*F 1(2, 1 ) + S* F 1 (2,2) 

F2 ( 1 , 2 ) = ~S* F 1 ( 1 , 1 ) - C*F 1 ( 1 , 2 ) 

F2 ( 2 , 2 ) = -S*F 1(2,1) - C*F1(2.2) 

RETURN 

END 

C 

SUBROUTINE FR I FORC ( DMATRX , SMATRX , FCOE # FNORM , FFMATRX ) 

C 

C CALCULATE FRICTION FORCE MATRIX WHEN RELATIVE 

C DISPLACEMENT MATRIX OF THE CONTACT POINTS 

C IS GIVEN 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION DMATRX ( 2 ,2), SMATRX ( 2 , 2 ) , FFMATRX ( 2 „ 2 ) 
DIMENSION SVALUE ( 2 ) , RANGLE ( 2 ) , PANGLE(2) ,U(2, 2) f V(2,2) 
DIMENSION FFSTIC(2 t 2) ,FFSLIP(2, 2) ,RFFSLIP(2,2) 

CALL AB(SMATRX, DMATRX, 2, FFSTIC) 

CALL ROT ATE(FF STIC, SVALUE , RANGLE , P ANGLE ) 

I F ( SVALUE ( 1 ) . EQ.0.0) THEN 
DO 10 1=1,2 

DO 10 J=1 ,2 

10 FFMATRX ( I , J) = 0.0 
GO TO 100 
ENDIF 

FLIMIT = FCOE* FNORM 
SLIP = FLIMIT/SVALUE( 1 ) 

I F ( SLI P . LT . 1.0) GO TO 20 
C 

CALL EMATRXCFFSTIC, FFMATRX, 2) 

C 

GO TO 100 
20 CONTINUE 
C 

CALL ROTATE( DMATRX , SVALUE , RANGLE , PANGLE ) 

R = SVALUE ( 2 ) / SVALUE ( 1 ) 

CALL ELLIPI (R.BSIN.BCOS) 

FOURDPI = 4.0/3.1415927 
RFFSLIPC 1,1) = 0.0 

RFFSL I P( 1 . 2 ) = -FLIMIT ♦ FOUR DP I *OS I N 
RFFSLIP(2, 1 ) = FLIMIT*FOURDPI *BCOS 

RFFSLI P( 2 , 2) = 0.0 
C 

CALL RTMATRX ( RANGLE ( 1 ) ,U) 

CALL RTMATRX ( -PANGLE ( 1 ), V) 

CALL ABC(U,RFFSLIP,V, 2 , FFSLIP) 

C 
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c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


INTERPOLATION COEFFICIENTS 

CSTHETA = 1.0 - 2 . 0 * S L I P 
THETA = ACOS(CSTHETA) 

CSTIC = (THETA - 0 . 5 * S I N ( 2 . 0* THETA ) ) / 3 . 1 4 1 59 27 
CSLIP = 1 . 0 - SLIP 

CALL COMATRX (CSTIC, FFSTIC , CSLIP, FFSLIP, 2 . FFMATRX) 
100 CONTINUE 
RETURN 
END 

FUNCTION ELL I P 1 (T) 


ELLIP1 - COMPLETE ELLIPTIC INTEGRAL OF THE FIRST KIND. 

T - INPUT PARAMETER; 0 < T < 1.0 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

AO = 1 .386294361 12 
A 1 = 0.09666344259 
A 2 = 0.03590092383 
A3 = 0.03742563713 
A4 = 0.01451 196212 
BO = 0.5 

B1 = 0. 12498593597 
B2 = 0.06880248576 
B3 = 0.03328355346 
B4 = 0.00441787012 
T1 = 1 .0 - T 
T 2 = T 1 * * 2 
T 3 = T 1 * *3 
T 4 = T 1 * *4 

ELL I P 1 = (A0 + A1 *T1+A2*T2 + A3*T3*-A4*T4) 

ELL I P 1 = ELL I P 1 + ( BO + B 1 *T1+B2*T2 + B3*T3 + B4*14 ) +LOG ( 1 . 0/ T 

RETURN 

END 

FUNCTION ELLIP2 (T) 

ELLIP2 - COMPLETE ELLIPTIC INTEGRAL OF THE SECOND KIND. 

T - INPUT PARAMETER; 0 < T < 1.0 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

AO =1.0 

A 1 = 0.44325141463 

A 2 = 0.06260601220 
A3 = 0.04757383546 
A4 = 0.01736506451 
B 1 = 0.24998368310 
B2 = 0.09200180037 
B3 = 0.04069697526 
B4 = 0.00526449639 
T1 = 1 .0 - T 
T2 = T 1 * * 2 
T3 = T 1 * * 3 
T4 = T 1 * * 4 

ELLIP2 = (AO+A 1 *T1+A2*T2+A3*T3+A4*T4) 

ELLIP2 = ELLIP2 + ( B 1 * T 1 + B 2 * T 2 + B3 * T 3 + B4 * T4 ) * LOG ( 1 . 0 / T 1 ) 

RETURN 


) 
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END 

C 

SUBROUTINE ROTATE ( A , SVALUE , RANGLE , PANGLE ) 

C 

C CALCULATE SINGULAR VALUES AND VECTORS OF A MATRIX A 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A (2, 2) ,SVALUE(2) , RANGLE ( 2 ) , PANGLE ( 2 ) 

C 


DS = 

A( 1 , 1) 

+ 

A ( 2 , 

,2) 

DD = 

A ( 1 , 1 ) 

- 

A ( 2 , 

, 2 ) 

CS = 

A( 1 ,2) 


A ( 2 , 

> 1 ) 

CD = 

A ( 2 , 1 ) 

- 

A ( 1 , 

2) 

SDS 

= DS*DS 




SDD 

= DD*DD 




SCS 

= cs*cs 




SCD 

= CD*CD 





SVALUE ( 1 ) = 0.5*( (SDS+SCD)**0.5 + ( SDD+ SCS ) **0.5) 

SVALUE ( 2 ) = 0 . 5* ( (SDS+SCD) **0.5 - ( SDD+SCS ) **0.5) 

I F ( SVALUE ( 1 ) . EQ . SVALUE ( 2 ) ) GO TO 10 
SRAP = ATAN2(CS,DD) 

GO TO 20 
10 SRAP = 0.0 

20 I F ( SVALUE ( 1 ) .EQ. (-SVALUE(2) ) } GO TO 30 
DRAP = AT AN2 ( CD ,DS) 

GO TO 40 
30 DRAP = 0.0 
40 CONTINUE 

RANGLE ( 1 ) = 0 . 5* ( SRAP + DRAP) 

RANGLE ( 2 ) = RANGLE ( 1 ) * 1 80 . 0 / 3 . 1 4 1 59 2 7 
PANGLE ( 1 ) = 0 . 5* ( SRAP - DRAP) 

PANGLE ( 2 ) = PANGLE( 1 )♦ 180.0/3. 1415927 

RETURN 
END 
C 

SUBROUTINE RTMATRX (RANGLE, U) 

C 

C GIVE TWO DIMENSIONAL ROTATION MATRIX WHEN ROTATIONAL 

C ANGLE IS GIVEN 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION U ( 2 , 2 ) 

U( 1 , 1 ) = COS (RANGLE ) 

U ( 2 , 1) = SIN (RANGLE ) 

U( 1 , 2) = -U( 2 , 1 ) 

U ( 2 , 2 ) = U( 1 , 1 ) 

RETURN 

END 

C 

SUBROUTINE ELL I P I ( R , BS I N , BCOS ) 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

I F ( ABS ( R ) .LE. 1.0D-04) GO TO 25 
IF(ABS(R) .EQ. 1.0) GO TO 20 
IF(ABS(R) .GT. 1.0) GO TO 10 
T = 1.0- R*R 

El = ELLIP1 (T) 

E2 = ELLIP2(T ) 

BSIN = ( E2 - (1 . 0-T ) *E 1 )/T 
BCOS = (El - E 2 ) *R / T 
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GO TO 30 
10 CONTINUE 

T = 1.0- 1 / ( R*R ) 

El = ELLIP1 (T) 

E 2 = ELLIP2(T) 

BCOS = (E2 - ( 1 . 0-T ) * E 1 )/T 
BSIN = (El - E2)/R/T 
GO TO 30 
20 CONTINUE 

BSIN = 3. 1415927/4.0 
BCOS = R * BS I N 
GO TO 30 
25 BSIN = 1.0 
BCOS = 0.0 
30 CONTINUE 
C 

RETURN 

END 

C 

SUBROUTINE AB(A,B,N,C) 

C 

C C = AB (MATRIX MULTIPLICATION) 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
DIMENSION A(2,2) , B ( 2 , 2 ) , C ( 2 , 2 ) 

DO 10 1=1 ,N 
DO 10 J = 1 , N 
C(I . J) = 0.0 
DO 20 K= 1 , N 

20 C(I,J) = C(I,J) + A ( I ,K)*B(K,J) 

10 CONTINUE 
C 

RETURN 

END 

C 

SUBROUTINE ABC ( A , B , C , N , D ) 

C 

C D = ABC (MATRIX MULTIPLICATION) 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A ( 2 . 2 ) , B ( 2 , 2 ) ,C(2,2) , D ( 2 , 2 ) ,£(2,2) 
CALL AB(A,B,N,E) 

CALL AB(E.C.N.D) 

RETURN 

END 

C 

SUBROUTINE EMATRX ( A , B , N ) 

C 

C DUPLICATES A MATRIX A = B 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A ( 2 , 2 ) ,B(2.2) 

DO 10 I = 1 . N 
DO 10 J=1 ,N 
10 B(I.J) = A ( I , J ) 

RETURN 

END 

C 

SUBROUTINE COMATRX ( C A , A , CB , B , N , C ) 

C 



C A + A + CB + B 



C LINEAR COMBINATION OF TWO MATRICES C = 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION A(2,2) f B(2,2) , C ( 2 , 2 ) 

DO 10 1 = 1 ,N 
DO 10 J = 1 , N 

10 C(I,J) = CA*A(I,J) + CB * B ( I ,J) 

RETURN 
END 
C 

SUBROUTINE GSLV ( NOEQN , X , DX , XMAX , XMI N , DELMX , EPS I , MITER, 

1 I PARTL , I ERR , A , B ) 

C 

C 

C GSOL V - SOLVES A NONLINEAR EQUATIONS SET BY NEW ! ON - R A PHSON ' S 

C METHOD. 

C 

C PROGRAM DESCRIPTION 

C SOLVES A NONLINEAR EQUATION SET BY NEW1 ON-RAPHSON ' S 

C METHOD. (UP TO 100 X 100 EQUATIONS SET). 

C 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COMMON/NITER/ITER 

DIMENSION X(NOEQN) , XMAX ( NOEQN ) , XMI N ( NOEQN ) , 

1 DELMX ( NOEQN ) , DX ( NOEQN ) , EPS I ( NOEQN ) 

DIMENSION A (NOEQN, NOEQN) , B ( NOEQN , 1 ) 

DIMENSION LW(300) ,XN( 100) ,FV( 10) 

EXTERNAL F 
C 

C INITIALIZE DATA 
C 

I ERR=0 
I HALT=0 
I TER=0 
C 

C SOLVE THE EQUATION USING THE NEWTON'S METHOD 
C 

10 CONTINUE 
C 

C IF NO. OF ITERATIONS ARE GREATER THAN MITER, DIVERGES 
C 

I TER= I TER+ 1 

IF ( ITER .GT.M ITER) THEN 
I ERR= 1 
GO TO 70 
ENDIF 
C 

C GET THE FUNCTION VALUE FOR THE CURRENT VARIABLE VALUES 
C 

DO 40 1=1 .NOEQN 
I ENUMB= I 

VALF= FOENUMB, NOEQN, X) 

8(1,1)= -VALF 
FV ( I ) =\/ALF 
C 

C GET THE PARTIAL DIFFERENCE VALUE FOR THE CURRENT VARIABLE VALUES 
C 

DO 30 J= 1 . NOEQN 
I VNUMB= J 
C 
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C STOP WHEN THE VALUE OF THE VARIABLE IS LESS THAN XMIN 
C OR GREATER THAN XMAX . 

C 

I F ( X ( I VNUMB ) . LT . XMI N ( I VNUMB ) ) THEN 
X ( I VNUMB )=X( IVNUMB)+500 . 

GO TO 10 

EL SEIF C X ( I VNUMB) . GT . XMAX( I VNUMB) ) THEN 
I ERR = 2 
GO TO 70 
ENDIF 
C 

C CALCULATE THE PARTIAL DERIVATIVES NUMERICALLY. 

C 

IF( IPARTL. EQ.O) THEN 
DO 20 K= 1 , NOEQN 
XN ( K ) = X ( K ) 

20 CONTINUE 

XN ( I VNUMB ) =X ( I VNUMB ) +DX ( I VNUMB ) 

FNEW = FOENUMB, NOEQN, XN) 

A ( I , J) = ( FNEW- VA LF ) /DX( 1VNUMB) 

C 

C CALCULATE THE PARTIAL DERIVATIVES ANALYTICALLY. 

C 

ELSEIF( IPARTL. EQ. 1 ) THEN 
C A(I,J)= PARTL ( I ENUMB , I VNUMB , N , X ) 

ENDIF 

C 

30 CONTINUE 
40 CONTINUE 
I HALT = 1 

DO 44 1=1 , NOEQN 
DAA=ABS(FV( I ) ) 

I F ( DA A . GT . EPSI ( I ) ) THEN 

I HALT = 0 

ENDIF 

44 CONTINUE 

I F ( IHALT. EQ. 1 ) GO TO 70 
C 

C SOLVE THE NOEQN BY NOEQN METRICS 
C 

CALL ULNEQ2( A, NOEQN, NOEQN. B, 1 .NOEQN, LW,0. . LRR ) 
IF(LRR.EQ.I) THEN 
I ERR=3 
GO TO 70 
ENDIF 
C 

C IF THE CHANGE OF VARIABLE IS GREATER THAN THE ALLOWABLE 
C MAXIMUM VALUE, LIMIT THE CHANGE OF THE VARIABLES. 

C 

DO 50 1=1 , NOEQN 
I VNUMB= I 

I F ( ABS ( B ( 1,1)) . GT . DELMX ( I VNUMB) ) THEN 
S I GN = B ( I . 1 ) /ABS ( B ( I . 1 ) ) 

B ( I , 1 ) =SI GN* DELMX ( I VNUMB) 

ENDIF 

50 CONTINUE 
C 

C GET THE NEW VALUES OF THE VARIABLES 
C 


IHALT= 1 
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I ERR = 2 
GO TO 70 
END I F 
C 

C CALCULATE THE PARTIAL DERIVATIVES NUMERICALLY. 

C 

I F ( I PARTL . EQ . 0 ) THEN 
DO 20 K= 1 , NOEQN 
XN ( K ) = X ( K ) 

20 CONTINUE 

XN ( I VNUMB ) =X ( I VNUMB ) +DX ( I VNUMB) 

FNEW = FI ( IENUMB, NOEQN ,XN) 

A ( I , J ) = ( FNEW-VALF ) /DX ( I VNUMB ) 

C 

C CALCULATE THE PARTIAL DERIVATIVES ANALYTICALLY. 

C 

ELSEIF(IPARTL,.EQ. 1 ) THEN 
C A ( I , J ) = PARTL ( I ENUMB , 1 VNUMB , N , X ) 

ENDIF 

C 

30 CONTINUE 
40 CONTINUE 
I HALT= 1 

DO 44 1=1 , NOEQN 
DAA=ABS ( F V ( I ) ) 

I F (DA A .GT . EPS I ( I ) ) THEN 

I HALT=0 

ENDIF 

44 CONTINUE 

I F ( I HALT . EQ . 1 ) GO TO 70 
C 

C SOLVE THE NOEQN BY NOEQN METRICS 
C 

CALL ULNEQ2(A , NOEQN , NOEQN, B, 1 , NOEQN, L.W.O. , LRR ) 
IF(LRR.EQ.I) THEN 
I ERR=3 
GO TO 70 
ENDIF 
C 

C IF THE CHANGE OF VARIABLE IS. GREATER THAN THE ALLOWABLE 
C MAXIMUM VALUE, LIMIT THE CHANGE OF THE VARIABLES. 

C 

DO 50 1=1 , NOEQN 
I VNUMB= I 

IF(A8S(B(I , 1 )) .GT. DELMX ( I VNUMB ) ) THEN 
SIGN=B( I , 1 ) / ABS ( B ( 1,1)) 

B ( I , 1 )=SIGN*DELMX( I VNUMB) 

ENDIF 

50 CONTINUE 
C 

C GET THE NEW VALUES OF THE VARIABLES 
C 

I HALT = 1 

DO 60 1=1 , NOEQN 
DELX = B ( I , 1 ) 

XNEW = X ( I ) +B ( I , 1) 

X ( I ) = XNEW 
C 

60 CONTINUE 
GO TO 10 
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70 RETURN 
END 
C 

SUBROUTINE ULNEQ2 ( A , N , I D I MA , B . M , I D I MB . LWORK , L PS , I ERR ) 
ULNEQ2 - SOLVES A SET OF DOUBL E - PR EC I S I ON LINEAR EQUATIONS. 


PROGRAM DESCRIPTION 

THIS SUBROUTINE SOLVES THE MATRIX EQUATION A * X = B , 

OVERWRITING B WITH THE SOLUTION MATRIX X. A MUST BE SQUARE 
AND NON-SINGULAR. B MUST HAVE THE SAME NUMBER OF ROWS AS A. 

BOTH A AND B ARE DESTROYED. BOTH A AND B ARE DOUBLE- 
PRECISION MATRICES. 

ACCESS 

CALL ULNEQ2 ( A . N , IDIMA , B ,M , IDIMB . LWORK . EPS , I ERR ) 

A -DP ARY-U/R-DIM( IDIMA, N)- COEFFICIENT MATRIX OF THE 

MATRIX EQUATION A * X = B. THE N X N PORTION OF A 
MUST CONTAIN THE SQUARE COEFFICIENT MATRIX. THE 
CONTENTS OF THIS PORTION ARE DESTROYED BY THIS ROUTINE. 

N -IN VBL-USD- THE ORDER OF COEFFICIENT MATRIX A . N MUST 

BE LESS THAN OR EQUAL TO IDIMA, LESS THAN OR EQUAL TO 
THE SECOND DIMENSION OF A IN THE CALLING ROUTINE, AND 
LESS THAN OR EQUAL TO IDIMB. 

IDIMA -IN VBL-USD- FIRST DIMENSION OF A IN THE CALLING ROUTINE. 

B -DP ARY-U/R-DIM( IDIMB ,M) - CONSTANT-TERM ( R I GHT-HAND- S I DE ) 

MATRIX. THE N X M PORTION OF B MUST CONTAIN THE 
RIGHT-HAND-SIDE MATRIX. UPON RETURN FROM THIS ROUTINE, 
THIS PORTION OF B CONTAINS THE SOLUTION MATRIX. 

M -IN VBL-USD- THE NUMBER OF COLUMNS IN B. M MUST BE LESS 

THAN OR EQUAL TO THE SECOND DIMENSION OF B IN THE 
CALLING ROUTINE. 

IDIMB -IN VBL-USD- FIRST DIMENSION OF B IN THE CALLING ROUTINE. 
LWORK -IN AR Y-WRK-D I M ( 3*N ) - WORKING ARRAY FOR SOLUTION 

ALGORITHM. LWORK MAY BE ANY INTEGER ARRAY WITH AT LEAST 
3N ELEMENTS. THE FIRST THIRD IS RESERVED FOR PIVOT 
COLUMN INDICES. THE FIRST NP POSITIONS OF THIS THIRD 
LIST THE PIVOT COLUMN INDICES IN ORDER OF USE. THE 
SECOND THIRD IS RESERVED FOR PIVOT ROW INDICES. THE 
FIRST NP POSITIONS OF THIS THIRD LIST THE PIVOT ROW 
INDICES IN ORDER OF USE. THE LAST THIRD IS USED FOR 
TEMPORARY STORAGE FOR INTERCHANGE OF P l VO i ROW AND 
COLUMN INDICES. 

EPS -RL VBL-USD- VALUE ALL PIVOT ELEMENTS MUST EXCEED FOR 

MATRIX A TO BE CONSIDERED NONSINGULAR. IF IN DOUBT, 
USE THE VALUE 0.0 . 

I err -IN VBL-RTD- PIVOT-SEARCH ERROR CODE (SEE BELOW). 

COMMON BLOCK VARIABLES 
NONE 

ERROR CONDITIONS 

PIVOT-SEARCH ERRORS ARE RETURNED THROUGH I ERR AS FOLLOWS- 

I ERR = 0 IF ALL COLUMNS OF X ARE FOUND , NO TROUBLE BEING DETECTED. 

IERR=1 IF NO COLUMNS OF X ARE FOUND, THE ELIMINATION PROCESS 

BEING HALTED BECAUSE THE CURRENT PIVOT FAILS TO EXCEED EPS IN 
MAGNITUDE. 

EXTERNAL REFERENCES 




to 

u> 

o 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


NONE 

COMMENTS 

THE METHOD CONSISTS OF GAUSSIAN ELIMINATION FOLLOWED BY BACK 
SUBSTITUTIONS. THIS IS MORE EFFICIENT THAN SOLUTION BY MATRIX 
INVERSION REGARDLESS OF THE NUMBER OF COLUMNS IN B. BOTH ROWS 
AND COLUMNS ARE SEARCHED FOR MAXIMAL PIVOTS. INTERCHANGING OF 
ROWS OR COLUMNS OF A IS AVOIDED. CHAPTER 1 OF E.L. STIEFLE, 
INTRODUCTION TO NUMERICAL MATHEMAT I C S , ACADEM I C PRESS , N . Y 1 963 , 
SHOULD BE HELPFUL IN FOLLOWING THE CODE. 

LOCAL VARIABLES 

ABSEPS “DP VBL- ABSOLUTE VALUE OF EPS. 

ABSPIV -DP VBL- ABSOLUTE VALUE OF PIV. 

IPIV -IN VBL- ACTUAL ROW OF CURRENT PIVOT ELEMENT. 

ITEMP -IN VBL- TEMPORARY SPACE FOR INTERCHANGE OF PIVOT 
INDICES. 

JPIV -IN VBL- ACTUAL COLUMN OF CURRENT PIVOT ELEMENT. 

KPIV -IN VBL- LOCATION IN SECOND THIRD OF LWOHK (ROW PIVOTS) 

CORRESPONDING TO CURRENT PIVOT. 

LPI V -IN VBL- LOCATION IN FIRST THIRD OF LWORK (COLUMN PIVOTS) 

CORRESPONDING TO CURRENT PIVOT. 

NP -IN VBL- NUMBER OF PIVOT ELEMENT CURRENTLY BEING 

COMPUTED. 

NPN -IN VBL- N PLUS N (I.E. 2N). 

NPP -IN VBL- NP PLUS 1 . 

PIV -DP VBL- VALUE OF THE CURRENT PIVOT. 

TEMP -DP VBL- TEMPORARY SPACE FOR INTERCHANGE OF ELEMENTS OF 
SOLUTION MATRIX. 

DIMENSION A ( I DIMA , N ) , B ( I DI MB , M ) .LWORK (300) 

DOUBLE PRECISION A , B 

DOUBLE PRECISION ABSEPS , ABSPI V , P I V , TEMP 


INITIALIZATIONS 

I ERR= 1 
NPN = N + N 

ABSEPS=ABS(EPS) 

DO 1 1 = 1 ,N 
LWORK ( I +N ) = I 
LWORK ( I )=I 
1 CONTINUE 

BEGIN ELIMINATION PROCESS 

DO 10 NP= 1 , N 

SELECT PIVOT 

ABSPI V=0 . 0 
DO 3 K=NP , N 
I =LWOR K ( K + N ) 

DO 2 L=NP , N 
J = LWOR K ( L ) 

IF (ABS(A(I ,J) ) .LE. ABSPIV) GOTO 2 
KPI V = K 
LP I V = L 
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IPIV=I 
JPIV=J 
PI V = A ( I , J ) 

ABSPIV=ABS(PIV) 

2 CONTINUE 

3 CONTINUE 

C 

C EXIT IF PIVOT TOO SMALL 

C 

C 

IF (ABSPIV.LE. ABSEPS) GOTO 19 
C 

C UPDATE PIVOT ROW AND COLUMN LISTS 

C 

I TEMP=LWORK ( NP+N ) 

LWORK ( NP+N ) = LWORK (KPIV+N) 

LWORK ( KPIV+N)=ITEMP 
I TEMP= LWORK ( NP ) 

LWORK ( NP ) =LWORK ( LPI V ) 

LWORK (LPIV)=ITEMP 

C MODIFY PIVOT ROW OF A AND B (ELEMENTS IN PRESENT OR PREVIOUS 

C PIVOT COLUMNS OF A ARE SKIPPED) 

C 

IF (NP.EQ.N) GOTO 5 
NPP=NP+ 1 
DO 4 L=NPP,N 
J = LWORK ( L ) 

A(IPIV, J)=-A( IPIV, J)/PIV 

4 CONTINUE 

5 DO 6 J= 1 ,M 

B( IPIV, J)=-B( IPIV t J) /PI V 

6 CONTINUE 
C 

C MOOIFY NON-PIVOT ROWS OF A AND B (ELEMENTS IN PRESENT OR 

C PREVIOUS PIVOT ROWS OR COLUMNS ARE SKIPPED) 

C 

IF (NP.EQ.N) GOTO 10 
DO 9 K=NPP , N 
I =LWORK ( K + N ) 

TEMP=A ( I , JPI V) 

IF (TEMP. EQ. 0.0) GOTO 9 
DO 7 L=NPP t N 
I -l WHRK ( 1 ! 

A CI .J)=A(I ,J)+A(IPIV.J)*TEMP 

7 CONTINUE 
DO 8 J=1 .M 

B(I,J)=B(I , J)+B(IPIV,J)*TEMP 

8 CONTINUE 

9 CONTINUE 
10 CONTINUE 

C 

C END ELIMINATION PROCESS 

C 

C BEGIN BACK SUBSTITUTIONS 

C 

DO 13 J=1 ,M 
DO 12 K = 2 , N 
KK=N-K+ 1 
I = LWOR K ( KK + N ) 
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DO 1 1 L-2 , K 

LL=N-L+2 
I I=LWORK(LL+N) 

JJ=LWORK(LL) 

B(I,J)=8(I , J ) + B ( I I , J ) * A ( I ,JJ) 

11 CONTINUE 

12 CONTINUE 

13 CONTINUE 
C 

C UNSCRAMBLE ROWS OF SOLUTION MATRIX 

C 

DO 14 1=1 ,N 
L=LWORK ( I +N ) 

LWORK(L + NPN) =LWOR K ( I ) 

14 CONTINUE 

DO 17 1=1 ,N 

15 K = LWOR K( I+NPN ) 

IF (I.EQ.K) GOTO 17 
DO 16 J=1 ,M 
TEMP=B ( I , J) 

B(I,J)=B(K,J) 

B ( K , J ) =TEMP 

16 CONTINUE 

LWORK( I+NPN )=LWORK( K+NPN) 

LWOR K(K + NPN)=K 
GO TO 15 

17 CONTINUE 
C 

DO 18 J= 1 , M 
DO 18 1=1 ,N 
8 ( I , J ) =-B ( I , J ) 

18 CONTINUE 
I ERR=0 

19 RETURN 
END 
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